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Reaction  bonded  silicon  nitride  (RBSN)  prepared  from  submicron  size  pure 
silicon  powder  is  a  micro-porous  ceramic  with  excellent  mechanical  and 
dielectric  properties  for  high  temperature  structural  applications  and  i 
for  radomes.  In  this  research  program  the  mechanical  properties  includ¬ 
ing  elastic  moduli,  tensile  strength  and  fracture  toughness  of  RBSN  was 
experimentally  studied  in  connection  with  its  quasi-regular  microstruct 
^ure.  The  further  ranges  of  microstructural  features  and  the  processes  . 
related  to  local  fracture  resistance  were  studied  by  means  of  two  relate 
and  quite  flexible  computer  simulations  of  the  reaction  kine'tics  control 
ing  the  microstructure  and  the  processes  that  control  the  crack  paths. 
These  results  are  summarized  and  presented  in  greater  detail  in  four 
appended  manuscripts.  \ 
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MECHANICAL  PROPERTIES  OF  POROUS  HIGH 
TEMPERATURE  STRUCTURAL  MATERIALS: 

SOURCES  OF  TOUGHNESS  IN  REACTION  BONDED  SILICON  NITRIDE 


1 . 0  Introduction 

In  the  production  of  high  temperature  ceramics  by  gaseous  reaction  bonding,  as  in  the  case 
of  reaction  bonded  silicon  nitride  (RBSN),  a  material  of  quasi-homogenous  porous 
microstructure  is  obtained  with  a  porosity  in  the  range  of  0.2  -  0.3.  A  ceramic  such  as 
RBSN,  prepared  with  no  liquid  phase  sintering  aids,  has  many  attractive  high  temperature 
properties  and  physical  characteristics.  It  is  a  candidate  material  for  the  structural 
components  of  the  Advanced  Turbine  Systems  with  potential  operating  temperatures  in  the 
range  of  1400C.  With  careful  processing  RBSN  can  be  produced  in  near  net  shapes.  In 
fully  reacted  form  with  a  minimum  of  retained  unreacted  Si  it  has  a  low  dielectric  constant 
and  exhibits  low  dielectric  losses  making  it  a  very  attractive  material  for  radome  applications. 
The  special  processing  procedures  of  this  material  and  its  structural  and  dialectic  properties 
that  have  been  developed  and  explored  under  previous  ARO  and  NASA  support  are 
summarized  in  Appendix  A. 


From  a  more  generic  point  of  view,  RBSN  of  porosity  in  the  range  of  0.25,  is  a  material 
nearly  an  order  of  magnitude  more  porous  than  the  levels  of  porosity  encountered  in 
imperfectly  sintered  solids,  but  is  very  much  denser  than  the  usual  cellular  solids.  While  it 
has  many  microstructural  characteristics  common  with  materials  on  both  sides  of  the 
porosity  scale,  RBSN  requires  different  methodology  for  understanding  of  its  mechanical 
properties.  In  the  present  case  it  was  concluded  that  the  attractive  structural  potential  and 
other  outstanding  physical  characteristics  of  RBSN  would  be  difficult  to  exploit  without  fully 
understanding  and  improving  its  mechanical  properties,  and  in  particular  its  strength  and 
fracture  toughness.  Thus,  in  the  AFOSR  program  while  considerable  effort  was  direction  to 
the  processing  aspects  of  RBSN  from  pure  specially  prepared  Si  powder,  the  main  emphasis 
was  on  the  understanding  of  its  mechanical  properties  and  particularly  the  factors  that  affect 
the  fracture  toughness.  In  this,  a  complete  cycle  approach  was  taken  as  a  case  study  in  the 
development  of  an  advanced  material  by  design.  In  the  complete  cycle  approach  to  an 
attractive  candidate  material  such  as  RBSN,  the  starting  point  is  a  detailed  study  of  the 
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physical  process  beginning  with  the  gaseous  reaction  processing,  to  the  topology  of  the 
reaction  induced  microstrueture  and  the  characteristics  of  this  microstructure  that  directly 
affect  its  strength  and  fracture  toughness.  This  is  followed  by  a  detailed  computer  simulation 
of  both  the  gaseous  reaction  process  controlling  the  microstructure  and  the  fracture  paths  that 
are  taken  by  a  crack  through  the  resulting  microstmcture.  When  satisfaetory  congmity  is 
established  between  the  actual  physical  phenomenon  of  the  observed  structure  property 
relation  and  that  in  the  computer  simulations,  then,  the  latter  is  used  to  explore  other  possible 
structure  property  relations  through  the  computer  simulation  to  deeide  on  new  avenues  for 
processing  to  improve  properties. 

2.0  The  Material  and  Its  Mechanical  Properties 

The  RBSN  produced  at  M.I.T.  is  made  from  a  pure  silane-produced  (SiH4)  supply  of 
perfectly  spherical  sub-micron  size  Si  particles  through  a  gaseous  reaction  process  with  N2 
resulting  in  Si3N4  with  a  quasi-homogeneously  dispersed  porosity  of  about  0.25.  The 

pores  have  a  size  distribution  ranging  from  roughly  70-700nm  with  aspect  ratios  in  the  range 
of  1. 0-2.0.  The  individual  grains  in  the  microstructure,  on  the  average,  have  a  size  roughly 
half  that  of  the  pores.  A  certain  fraction  of  the  pores  are  intereonnected.  Uncharacteristically 
large  pores  of  large  aspect  ratio  which  could  strongly  reduce  the  strength  are  rare.  There  is 
x-ray  line  broadening  evidence  that  there  are  short  range  reaction  induced  residual  misfit 
stresses  in  the  microstructure  of  a  wave  length  roughly  equal  to  the  grain  size.  While  the 
specific  effect  of  these  residual  stresses  has  not  been  determined,  they  are  likely  to  be  a 
principal  cause  of  the  predominantly  intergranular  nature  of  the  fracture  in  the  RBSN. 

The  elastic  properties  of  all  RSBN  samples  of  various  porosities  were  measured  by  a  laser 
pulse-echo  technique  with  wave  length  much  larger  than  the  microstructural  seale.  The 
measured  Young’s  moduli  were  predictably  related  to  the  porosity.  The  typical  material  has 
strength  values  in  the  range  of  400MPa,  with  maximum  values  reaching  as  high  as  870MPa. 
Fracture  toughness,  ranged  from  2.2-2. 8MPa  Vm.  These  fracture  toughnesses  have  been 
obtained  from  fracture  loads  and  estimates  of  fracture  mirrors.  In  the  course  of  the  study  a 
new  method  has  been  devised  which  is  based  on  direct  measurement  of  shapes  of  open 
cracks  on  the  verge  of  propagation.  These  measurements  which  are  considered  more  reliable 


have  resulted  in  somewhat  lower  values  of  fracture  toughness.  In  such  determinations  of 
crack  shapes  some  evidence  of  crack  bridging  by  native  heterogeneities  in  the  microstructure 
has  been  found  in  addition  to  a  direct  correlation  between  crack  plane  tortuosity  and  level  of 
fracture  toughness.  These  measurements  and  features  were  taken  as  principal  characteristics 
of  the  RBSN  to  be  matched  by  simulations.  The  micro-structure  of  the  MIT  RBSN,  its 
stereological  description,  and  associated  mechanical  properties  have  been  described  in  detail 
by  Haubensak  et  al,  (1995a)  while  the  new  method  of  fracture  toughness  determination  has 
been  described  by  Haubensak  and  Argon  (1995).  These  manuscripts  are  attached  as 
Appendixes  B  and  C. 

3.0  The  Reaction  Kinetics  Simulation 

The  first  phase  of  the  simulation  involved  devising  a  reaction  kinetics  model  starting  with 
spherical  single  crystal  Si  particles  which  become  transformed  into  RBSN  through  a  gaseous 
reaction  with  N2.  The  reaction  kinetics  model  has  followed  a  cellular  automata  approach  in 
2-D,  with  some  crucial  topological  input  characteristic  of  3-D  features  of  the  actual  process. 

The  model  starts  by  considering  dense  randomly  packed  spherical  particles  of  constant  size 
in  space.  The  actual  coordinates  of  particle  centers  and  particle  sizes  for  the  2-D  model  are 
obtained  from  a  planar  cut  through  this  3-D  arrangement.  The  material  fraction  in  the  planar 
cut  is  then  transcribed  onto  a  hexagonally  tessellated  2-D  plane  with  mesh  size  roughly  0. 1 
that  of  the  disk  shaped  particles.  The  surface  sites  of  the  circular  particles  are  then  assumed 
to  fully  communicate  with  an  infinite  supply  of  N2  and  a  stochastic  gas  phase  reaction  with 
an  average  rate  proportional  to  the  Si  vapor  pressure  in  the  cavities.  The  overall  conversion 
reaction  is  governed  by  the  initial  density  of  Si3N4  nuclei  on  the  surfaces,  on  the  ratio  of 
reaction  rate  to  Si  evaporation  rate,  and  on  the  ratio  of  the  Si  evaporation  rate  to  the  re¬ 
condensation  rate.  Through  the  manipulation  of  these  parameters  and  the  incorporation  of 
some  surface  smoothing  by  diffusion  to  control  pore  surface  roughness,  computer  generated 
microstructures  with  features  closely  resembling  those  of  the  real  microstructures  were 
obtained.  In  this  simulation  a  rare  but  regularly  encountered  feature  was  some 
uncharacteristically  large,  ill  shaped  pores,  resulting  from  initial  packing  imperfections. 
These  are  expected  to  control  the  strength  and  are  considered  as  an  anologous  feature  to  the 


elongated  crack  like  imperfections  found  in  some  experimentally  produced  microstructures 
described  by  Haubensak,  et  al  (1995a).  On  the  whole,  the  reaction  kinetics  simulation  was 
successful  in  providing  realistic  rate  constants  for  the  overall  reaction,  for,  the  distribution  of 
grain  sizes  and  for  the  pore  size  and  their  aspect  ratio  distributions.  An  important 
observation  of  potential  relevance  for  improvement  of  mechanical  properties  was  the 
development  of  some  interpenetrating  elongated  grains  resulting  from  choosing  anisotropic 
reaction  directions.  These  were  observed  to  result  in  stronger  and  tougher  microstructures. 
This  simulation  has  been  described  in  detail  by  Haubensak  et  al  (1995b)  and  is  attached  to 
this  report  as  Appendix  D. 

4.0  The  Fracture  Path  Simulation 


The  fracture  path  simulation  was  also  performed  in  2-D  where  the  previously  simulated 
microstructure  is  overlaid  by  a  hexagonal  network  of  nodal  points  connected  by  linear  central 
force  springs  and  angle  flexing,  so  called,  “watch  springs”  with  spring  constants  chosen  to 
obtain  an  exact  match  to  the  elastic  properties  of  the  isotropic  RBSN.  The  nodes  and  the 
bond  springs  that  fall  into  pore  regions  of  the  microstructure  were  assigned  zero  strength  and 
zero  stiffness,  other  bonds  falling  into  the  grain  interiors  of  the  solid  matrix  material  were 
assigned  a  constant  common  stiffness,  while  those  bonds  traversing  across  grain  boundaries 
were  assigned  only  0.75  of  the  strength  of  matrix  bonds  to  mimic  the  intergranular  fracture 
behavior  of  the  real  material.  Quantitative  evaluations  were  based  on  a  bond-node-level 
stress  tensor  which  has  become  widely  used  in  simulations  of  heterogeneous  solids. 

Several  levels  of  strength  and  toughness  simulations  were  carried  out  with  the  model.  In  the 
first  level,  relatively  large  domains  of  typical  porous  microstructure  were  subjected  to  overall 
uniaxial  tensile  extension.  At  each  level  of  externally  imposed  displacements,  internal 
equilibrium  was  achieved  by  a  conjugate  gradient  method  of  energy  minimization.  Typical 
loading  sequences  consisted  of  searching  for  the  most  highly  stressed  bond,  releasing  it, 
incrementing  the  border  displacements  to  fracture  the  next  highest  stressed  bond  in  the  re¬ 
equilibrated  structure,  and  so  on,  until  the  entire  structure  could  no  longer  support  load. 
From  such  simulations  correlations  between  stiffness  and  peak  strength  of  representative 
volume  elements  (RVE)  were  obtained. 


In  a  higher  level  of  simulation  probing  the  crack  growth  resistance  of  such  microstructures 
the  entire  field  was  homogenized  in  density  but  each  hexagonal  node  and  its  environment 
were  assigned  a  randomly  selected  pair  of  strength  and  stiffness  values  obtained  from  the 
above  cell  simulations  of  RVEs.  Such  continuum  fields  with  random  properties  were  then 
probed  by  cracks  subjected  to  external  Mode  I  fields,  to  note  the  distribution  of  local  fracture 
toughness  and  the  resulting  tortusoity  of  the  crack  plane.  Again,  a  reasonably  close 
correspondence  was  obtained  between  the  actual  fracture  toughness  and  crack  plane 
tortusotities,  and  those  simulated  in  the  model.  Various  forms  of  microstructural  variability 
were  explored  in  the  simulation.  While  interesting  improvements  for  fracmre  toughness  was 
found  from  uncorrelated  increases  in  the  variability  of  bond  strength  (with  constant  bond 
stiffness)  as  well  as  with  bond  stiffness  (with  constant  bond  strength),  however,  for  the 
realistic  cases  where  the  strength  and  stiffness  were  correlated,  increasing  variability  resulted 
in  systematic  degradation  of  toughness. 

Finally,  3-D  features  in  the  local  crack  advance  were  approximated  by  the  introduction  of 
special  elements  that  incorporated  the  average  properties  of  a  random  sampling  of  three 
adjacent  2-D  elements,  as  a  first  approximation  of  a  typical  element  of  a  3-D  crack  front. 
Simulations  performed  with  such  elements  showed  more  realistic  variations  in  local  fracture 
toughness. 

The  entire  set  of  simulations  of  fracture  paths  have  been  described  by  Haubensak  et  al 
(1995c)  and  is  attached  to  this  report  as  Appendix  E. 

5.0  Conclusions 

From  such  simulations  of  the  microstructure  and  its  relation  to  strength  and  fracture 
toughness,  it  was  concluded  that  for  microstructure,  with  quasi-random  pore  homogeneity, 
substantial  improvement  of  toughness  is  not  likely.  Variation  of  grain  aspect  ratio  and 
interpenetration  of  elongated  grains  did,  however,  result  in  considerable  increases  in  fracture 
toughness.  This  suggests  that  attractive  improvements  in  the  fracture  toughness  of  such 
porous  solids  like  RBSN  may  still  be  possible  through  the  incorporation  of  a  randomly 
dispersed  population  of  whiskers. 


In  any  case  the  complete  cycle  approach  to  development  of  advanced  materials  by  design 

where  the  coupling  of  the  processing  with  property  simulations  is  accomplished  should  give 

superior  results. 

6.0  Personnel  Associated  with  the  Research  Program 

1 .  A.S.  Argon,  Principal  Investigator,  Quentin  Berg  Professor  of  Mechanical 
Engineering. 

2.  J.S.  Haggerty,  Co-Principal  Investigator,  Senior  Research  Scientist,  Materials 
Processing  Center. 

3 .  V.V.  Bulatov,  Research  Associate,  Department  of  Mechanical  Engineering  (part 
time). 

4 .  F.G.  Haubensak,  Research  Assistant,  Department  of  Materials  Science  and 
Engineering. 

5 .  A.  Lightfoot,  Staff  Researcher,  Materials  Processing  Center  (part  time). 

6 .  A.  Seman,  U.  Callinan,  Secretaries,  Mechanieal  Engineering  Department. 

7.0  Technical  Presentations 

1 .  A.S.  Argon,  “Interface  Fracture  and  Toughness  of  Composites”  Stanford  University 
Seminar  Lecture,  Nov.  19,  1992. 

2.  A.S.  Argon,  “Controlling  Toughness  in  Brittle  Composites”,  Lawrence  Livermore 
Laboratory  Lecture,  Nov.  6,  1992. 

3 .  A.S.  Argon,  “Controlling  Toughness  in  Brittle  Composites”,  Invited  lecture  at  the 
Institute  of  Mechanics  and  Materials  Symposium  entitled  “Interfaces  in  Materials  - 
Bridging  the  Length  Scale  Gap”,  La  Jolla,  CA,  Sep.  12,  1993. 

4.  A.S.  Argon,  V.V.  Bulatov,  and  F.G.  Haubensak,  “Simulation  of  the  Microstructure 
and  Fracture  Paths  of  RBSN”,  Contributed  presentation  at  the  Pacific  Rim  Meeting 
of  the  American  Ceramics  Society,  Honolulu,  Hawaii,  Nov.  8,  1993. 

5 .  A.S.  Argon,  “Toughening  Mechanisms  in  Brittle  Solids”,  Materials  Science  Seminar 
Lecture  at  Cornell  University,  Nov.  17,  1993. 

6 .  A.S.  Argon,  “Simulation  of  Crack  Propagation  Resistance  in  Brittle  Media  of  High 
Porosity”,  Invited  lecture  at  the  AIME  Symposium  on  High  Temperature  Materials  in 
Chicago,  IL,  Oct.  3,  1994. 


7 .  A.S.  Argon,  “Simulation  of  Crack  Propagation  Resistance  in  Brittle  Media  of  High 
Porosity”,  Seminar  lecture  at  the  Colorado  School  of  Mines,  Nov.  18, 1994. 

8 .  A.S.  Argon,  “Sources  of  Toughness  in  Reaction  Bonded  Silicon  Nitride”,  Applied 
Mechanics  Seminar  lecture.  Brown  University,  May  2,  1995. 

8.0  Theses  Completed 

Haubensak,  F.G.,  “Microstructure  Design  of  Porous  Brittle  Materials”,  Ph.D. 
Thesis,  Department  of  Materials  Science  and  Engineering,  MIT  September,  1994. 

9.0  Publications 

1 .  Haubensak,  F.G.  ,  Lightfoot,  A.,  Argon,  A.S.,  and  Haggerty,  J.S.,  (1995a) 
“Reaction  Bonded  Porous  Pure  Silicon  Nitride:  Microstructure,  Tensile  Strength, 
Fracture  Toughness”,  J.  Materials  Science,  submitted  for  publication. 

2.  Haubensak,  F.G.,  and  Argon,  A.S.,  (1995)  “A  New  Method  of  Fracture  Toughness 
Determination  in  Brittle  Ceramics  by  Open  Crack  Shape  Analysis”,  J.  Materials 
Science,  submitted  for  publication. 

3.  Haubensak,  F.G.,  Bulatov,  V.V.,  and  Argon.  A.S.,  (1995b)  “A  Simulation  of  the 
Reaction  Bonding  of  Si3N4  to  Generate  Porous  Microstructures”,  J.  Computer 
Aided  Materials  Design,  submitted  for  publication. 

4.  Haubensak,  F.G.,  Bulatov,  V.V.  and  Argon,  A.S.,  (1995c)  “Simulation  of  Crack 
Propagation  Resistance  in  Brittle  Media  of  High  Porosity”,  J.  Computer  Aided 
Materials  Design,  submitted  for  publication. 


Appendix  A;  Summary  of  Processing  Details  and  Early  Results  on 
Dielectric  and  Mechanical  Properties  of  Reaction  Bonded  Si3N4. 

(work  performed  under  earlier  ARO  and  NASA  support) 

A  two  year  program  designed  to  develop  high  purity  reaction  bonded  silicon  nitride  (RBSN) 
for  hypersonic  radome  applications  has  been  conducted  by  Raytheon  Corporation  and  the 
Massachusetts  Institute  of  Technology.  The  planned  three-year  program  was  terminated  one 
year  prematurely  by  funding  restrictions. 

The  development  program  was  initiated  based  oh  properties  exhibited  by  laboratory  scale 
RBSN  samples  made  by  MIT  from  semiconductor  grade  silane  (SiH4).  At  the  program’s 

inception  extraordinary  strengths  had  been  achieved  (averages  up  to  592  MPa  and  individual 
samples  up  to  870  MPa)  for  a  brittle  material  containing  ~  25%  porosity.  High  temperature 
dielectric  measurements  of  these  RBSN  samples  at  35  GHz  showed  low  losses  up  to 

1200®C  (tan  5  <  10'2),  low  dielectric  constants  (~  5.0)  and  low  temperature  coefficients. 

Net  shape  processing  experiments  demonstrated  that  nitrided  parts  exhibited  dimensions 
within  0.1 -0.2%  of  die  dimensions. 

The  goal  of  this  program  was  to  develop  a  process  for  fabricating  engineering  shapes 
(radomes)  from  this  material  which  retained  the  advantageous  properties  demonstrated  with 
laboratory  scale  samples,  and  to  acquire  a  data-base  of  mechanical,  dielectric  and  thermal 
properties  necessary  for  utilization  of  this  material.  Resulting  properties  were  analyzed  to 
determine  their  effect  on  permissible  flight  profiles. 

Alternative  forming  techniques  were  evaluated  in  terms  of  satisfying  the  dimensional  and 
quality  requirements,  and  also  in  terms  of  matching  existing  industrial  experience  and 
equipment  for  forming  large  radomes.  Considered  processes  included  dry  pressing,  isostatic 
pressing,  pressure  slip  casting,  colloidal  pressing,  centrifugal  casting,  electrophoretic 
deposition,  injection  molding  and  bulk  freeze  injection  molding.  Laboratory  scale 
experiments  demonstrated  the  feasibility  of  all  of  these  techniques  (except  for  bulk  freeze 
injection  molding).  Extensive  nitriding  kinetics  experiments  were  conducted  in  parallel  to 
insure  that  candidate  binders,  solvents  and  exposures  did  not  inhibit  complete  conversion  of 
the  Si  to  Si3N4  or  leave  residues,  both  of  which  are  essential  for  meeting  dielectic  property 


requirements.  The  combination  of  xylene  as  a  solvent  and  polystyrene  as  a  binder  proved 
satisfactory  in  terms  of  nitriding  rates,  extent  of  conversion,  residues  and  both  safety  and 
emission  criteria.  Powders  dispersed  in  this  solvent  and  binder  combination  were 
satisfactorily  spray  dried  to  form  handleable  powders.  Based  on  a  practical  solvent  and 
binder  combination,  a  demonstrated  powder  pretreatment  technique,  prior  radome  forming 
experience,  and  existing  industrial  forming  equipment,  warm  isopressing  of  spray  dried 
powders  was  selected  as  the  forming  process. 

Processing  experiments  were  initiated  at  MIT.  When  satisfactory  procedures  were 
identified,  duplicate  samples  were  made  at  Raytheon.  Extensive  characterizations  were  used 
to  demonstrate  the  equivalency  of  parts  made  in  both  laboratories.  Mechanical  and  dielectric 
properties  of  laboratory  scale  samples  made  by  warm  isopressing  were  made  to  qualify  the 
process  and  to  begin  optimization  of  process  variables. 

The  strengths  of  the  initial  batch  of  samples  made  by  warm  forming  spray  dried  Si  powders 
averaged  -308  MPa.  While  less  than  the  maximum  values  exhibited  by  optimized  colloidally 
pressed  samples  based  on  alcohol  solvents,  these  strengths  were  significantly  higher  than  is 
exhibited  by  conventional  RBSN  (-100-200  MPa)  and  much  higher  than  is  exhibited  by 
fused  silica  radome  material  (-48  MPa).  Later  generations  of  RBSN  samples  made  by  warm 
isopressing  spray  dried  powders  exhibited  slightly  improved  strengths  (335  MPa).  Using 
the  same  binder  and  solvent,  colloidally  pressed  samples  exhibited  average  strengths  of  539 
MPa.  These  strength  results  for  the  two  forming  techniques,  fractography  and 
characterization  of  the  spray  dried  powders  indicate  that  strengths  probably  were  controlled 
by  relics  of  oversize  particles  that  formed  during  spray  drying.  This  program  demonstrated 
that  high  strengths  are  achievable  with  the  selected  solvent  and  binder  through  uniformity  of 
microstructure.  Also,  strength  controlling  defects  can  be  removed  by  modifying  the  spray 
drying  process  and/or  removing  the  oversize  particles. 

Initial  high  temperature  dielectric  measurements  of  several  batches  of  warm  isopressed 
RBSN  by  the  Rockwell  Science  Center  gave  completely  different  results  than  the  colloidally 
pressed  sample  they  measured  prior  to  the  inception  of  the  program.  Resolving  the  cause  of 
these  anomalous  results  required  assembly  of  an  apparatus  to  make  elevated  temperature  (T  < 
500OC)  loss  tangent  measurements  at  Raytheon,  many  dielectric  measurements  and  many 


chemical  analyses.  Ultimately  it  was  found  that  the  high  losses  were  caused  by 
contamination  from  cements  used  to  mount  the  samples  at  Rockwell.  Subsequent 
measurements  showed  that  all  explored  process  conditions  produced  RBSN  samples  with 

loss  tangents  <5  x  lO"^  for  temperatures  up  to  SOO^C.  Also  important  for  radome 
applications,  the  dielectric  constants  were  low  (5.0-5.5)  and  relatively  insensitive  to 
temperature.  This  phase  of  the  program  demonstrated  that  improved  dielectric  properties 
were  reliably  achieved  for  RBSN  made  with  the  selected  procedures.  Although  optimization 
and  exploration  of  higher  temperature  properties  could  not  be  completed  within  the  shortened 

program,  all  of  the  data  indicates  that  acceptable  loss  tangents  (<  lO"^)  can  be  retained  for 
temperatures  in  excess  of  1200OC. 

Prototype  scale  samples  were  made  to  demonstrate  the  feasibility  of  making  practical  parts. 
Before  making  large  parts,  extensive  nitriding  kinetics  experiments  were  undertaken  to 
define  firing  schedules  which  avoided  excess  exothermic  heating  in  large  scale  parts  and 
which  also  yielded  residual  Si  levels  less  than  0.1%.  Using  spray  dried  powders,  and 
forming  conditions  defined  with  laboratory  scale  samples,  two  12-inch  long,  4-inch 
diameter,  3/16  inch  wall  cones  were  pressed  and  nitrided.  Both  of  the  as-formed  green  parts 
were  crack  free,  dimensionally  true  and  sound.  After  nitriding,  both  parts  exhibited  cracks 
and  evidence  of  overheating.  Irregularities  were  experienced  with  both  nitriding  runs  since 
these  represented  the  first  two  runs  in  a  much  larger  scale  furnace  than  had  been  used  with  all 
of  the  laboratory  scale  samples.  Detailed  comparisons  indicated  that  differences  in  firing 
schedules  and  differences  in  atmospheres  were  responsible  for  the  cracking.  No  evidence 
suggested  that  the  problem  was  intrinsic. 

Stress  levels  for  a  Patriot  configured  radome  subjected  to  a  PAC-3  flight  profile  were 
analyzed  by  Raytheon  Missile  Systems  Division.  Thermally  and  pressure  induced  stresses  at 
the  tip  and  the  base  of  the  radomes  were  determined.  Material  properties  were  based  on 
measured  and  projected  values  for  high  purity  RBSN  made  with  the  selected  binders, 
solvents  and  nitriding  schedules.  It  was  found  that  the  reduced  modulus  and  relatively  high 
thermal  conductivity  of  the  high  purity  RBSN  resulted  in  reduced  thermally-induced  stress 
levels  compared  with  fully  dense  Si3N4.  Total  induced  stress  levels  were  acceptable  in 

terms  of  measured  strengths. 


Although,  shortened,  this  program  accomplished  critical  objectives.  Forming  teehniques, 
binders,  solvents,  and  nitriding  schedules  were  identified  which  permitted  large  scale  RBSN 
parts  to  be  formed  from  high  purity  Si  powders.  The  property  levels  exhibited  by  laboratory 
scale  samples  made  with  the  selected  procedures  exhibited  properties  that  essentially  equaled 
the  best  levels  achieved  with  processes  designed  to  maximize  properties  without  regard  to 
practical  considerations.  These  outstanding  property  levels  will  permit  flight  profiles  that 
induce  significantly  higher  stagnation  temperatures  than  is  permitted  by  other  known 
materials.  The  elimination  of  machining  steps  should  also  result  in  much  reduced  production 
costs  -  not  to  mention  avoiding  machinery  induced  surface  damage.  Development  to  full 
scale  appears  practical  and  advantageous. 
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Abstract 


Porous  silicon  nitride  from  high  purity  uniform  sized  spherical  particles  of  sub-micron 
dimensions  obtained  by  a  gaseous  reaction  and  resulting  in  material  in  an  intermediate 
porosity  regime  of  0.2  -  0.3  was  studied  in  considerable  detail  to  establish  connections 
between  the  porous  microstructure  and  the  principal  mechanical  properties  such  as  Young’s 
Modulus,  tensile  strength,  and  fracture  toughness. 

In  the  course  of  the  investigation  several  new  techniques  of  material  characterization  were 
developed  which  proved  to  be  superior  to  what  has  been  used  in  the  past.  These  include  a 
laser  pulse-echo  technique  of  determining  the  Young’s  Modulus,  and  a  new  method  for 
determination  of  fracture  toughness  which  relies  on  accurate  measurements  of  the  shape  of 
cracks  at  impending  states  of  propagation. 

The  present  study  and  a  pair  of  related  computer  simulations  of  the  principal  elements  of  the 
gaseous  reaction  that  results  in  the  porous  microstructure  and  the  modes  of  the  crack 
propagation  through  such  microstructures  indicates  that  such  materials  with  rather  uniform 
porosity  have  a  limited  potential  for  development  of  tough  material.  Higher  fracture 
toughness  requires  the  introduction  of  new  degrees  of  freedom  that  nurture  crack 
deflection,  development  of  associated  redundant  cracking,  and  crack  bridging. 


*  Author  for  correspondence. 


1.0  INTRODUCTION 


The  cause-effect  relationship  between  microstructure  and  mechanical  properties  of  materials 
has  been  of  interest  for  several  decades.  In  particular,  the  characterization  of  the 
mechanical  behavior  of  heterogeneous,  and  specifically,  porous  materials,  has  recently 
become  of  technological  interest.  Research  has  typically  concentrated  on  materials  with 
highly  porous  cellular  microstructures  [1]  with  relative  densities  as  low  as  0.3%  and  on 
materials  with  dilute  porosity.  For  each  extreme  of  porosity  there  are  models  of  mechanical 
behavior  with  good  predictive  properties  that  are  based  on  engineering  assumptions  about 
the  microstructural  morphology.  In  each  case,  extrapolations  can  be  made  to  intermediate 
porosities  with  some  limited  success,  however,  there  is  a  fundamental  lack  of 
understanding  of  behavior  of  materials  within  the  intermediate  range  of  porosities. 

In  the  cases  where  the  effects  of  porosity  have  been  investigated,  often  the  study  relates  to 
relatively  isolated  pores  (and  consequently  low  levels  of  porosity)  acting  as  strength 
reducing  defects.  On  the  opposite  end  of  the  porosity  spectrum,  Gibson  and  Ashby  [1] 
have  assessed  the  field  of  cellular  materials  (p  <  30%).  These  density  trends  in  properties 
can  be  extended  somewhat  beyond  their  applicable  range  with  reasonable  results.  In  the 
intermediate  porosity  range,  the  mechanical  properties  of  ceramics  are  often  described  in 
terms  of  an  empirical  exponential  function  of  porosity.  Figure  1  schematically  illustrates 
the  limits  of  intermediate  porosity.  Rice  [2]  has  suggested  that  the  fracture  toughness  and 
moduli  of  low  porosity,  brittle  materials  can  be  described  by  such  an  exponential 
formulation.  While,  these  descriptive  equations  can  be  fitted  to  experimental  data,  they  do 
not  provide  pertinent  information  for  the  improvement  of  properties  through  microstructural 
control. 

Reaction  bonded  silicon  nitride,  developed  specifically  in  the  1950s  and  1960s  in  the  search 
for  high  strength  high  temperature  materials  [3],  is  intrinsically  porous.  Loosely  packed 
fine  silicon  particles  heated  at  temperatures  in  the  region  of  1 300C  in  a  stream  of  flowing 
nitrogen  were  found  to  convert  into  porous  silicon  nitride  with  considerable  strength.  The 
significant  processing  economies  of  the  low  processing  temperature  of  the  reaction  bonding 
process  relative  to  standard  sintering  processes  are  apparent  [4,5]  when  contrasted  with 
processes  that  produce  fully  dense  silicon  nitride,  such  as  including  hot  hydrostatic 
compaction  in  the  processing  steps  or  adding  densification  additives.  The  above 


processing  temperature  is  often  lower  than  the  target  service  operating  temperature  for  high 
temperature  diesel  or  gas-turbine  engine  hot  zone  components \  Nevertheless,  these 
beneficial  reaction  kinetics  allow  for  the  incoiporation  of  SiC  whiskers  with  minimal 
interaction  with  the  matrix  and  permit  complete  nitriding  in  less  than  10  hours  at  1400C  [6]. 
In  addition,  in  dynamic  applications  such  as  turbine  blades,  where  the  applied  stress  is 
linearly  proportional  to  the  mass  of  the  material  in  motion,  the  inherent  low  density  of 
porous  materials  translate  into  lower  operating  stresses.  Another  practical  advantage  is  the 
fact  that  many  reaction  bonding  processes  have  nearly  negligible  dimensional  change 
during  the  reaction  phase  due  to  the  fact  that  little  or  no  densification  by  sintering  occurs. 
Consequently,  cracking  associated  with  differential  shrinkage  is  avoided,  and  the  produced 
part  will  have  very  nearly  the  same  size  and  shape  as  the  precursor  green  body.  This  is  an 
advantage  if  the  produced  part  can  be  spared  the  final  machining  steps,  which  often 
introduce  surface  defects  that  are  well  known  to  weaken  ceramic  parts. 

However,  as  a  price  for  the  above  benefits,  the  final  microstructure  is  left  with  a  residual 
porosity.  The  consequence  of  this  porosity,  as  a  first  order  effect,  is  to  reduce  the  modulus 
and  fracture  energy,  both  of  which  have  been  observed  to  be  directly  related  to  the  volume 
fraction  of  porosity:  the  modulus  goes  to  zero  for  an  ideally  random  distribution  of 
porosity  at  approximately  15%,  the  three  dimensional  elastic  percolation  limit^  and  the 
fracture  energy  is  linearly  proportional  to  the  fracture  surface  area,  which  for  a  planar 
fracture  surface  implies  a  linear  dependence  on  porosity.  Thus,  while  mechanical 
properties  are  usually  inferior  at  high  porosities,  materials  have  been  produced  that,  on 
mass  specific  bases,  have  fracture  properties  that  are  comparable  to  fully  dense  materials 
[4].  There  are  some  indications  that  the  fracture  energy  could  be  significantly  increased  at 
higher  porosities  due  to  induced  crack  branching  [7]  and  grain  bridging  as  found  in  fully 
dense,  heterogeneous  material  [8,9].  The  RSBN  produced  at  MIT  with  high  purity  sub¬ 
micron  (0.3p,m)  diameter  Si  particles  from  laser-heated  SiH4  have  shown  unusually  rapid 

nitriding  rates  at  low  temperatures,  and  the  resulting  strengths  approach  those  of  fully 
dense  material  [5]. 


^  Gas  inlet  temperatures  approach  1370C  for  engine  designs  suitable  for  trucks  and  automobiles[3] 

*  The  critical  volume  fraction  at  the  percolation  threshold  for  a  variety  of  three  dimensional  lattices  ranges 
from  0.143  to  0.163.  Although  the  critical  volume  fraction  in  two  dimensions  based  on  continuum 
percolation  considerations  (exactly  0.5)  cannot  be  extended  directly  to  three  dimensions,  the  conjecture  that 
the  critical  volume  fraction  is  comparable  to  the  three  dimensional  packing  of  spheres  (=0. 1 5)  is  consistent 
with  Monte  Carlo  computation  (see  Ziman,  [10]). 


Previous  methods  to  determine  small  fracture  toughness  differences  of  RBSN  in  laboratory 
scale  specimens,  in  order  to  better  provide  for  material  optimization,  have  proved  to  be  less 
than  conclusive. 

Small  variations  in  fracture  toughness  by  standard  methods  were  not  detectable.  Ball-on- 
ring  strength  testing  and  corresponding  fractographic  determination  of  the  associated 
strength  limiting  flaw  geometry  proved  to  have  some  accuracy  but  was  affected  by  the 
subjective  nature  of  the  measurements.  While  they  serve  as  ready  indicators  of  broad 
trends  standard  indentation  fracture  toughness  techniques,  do  not  have  strong  fundamental 
accuracy,  and  uncertainties  associated  with  the  sensitivity  of  the  toughness  measurement  on 
the  crack  length  dimension  were  found  to  be  unacceptably  large.  To  rectify  these 
deficiencies,  a  technique  based  on  stable  crack  growth  from  an  internally  pressurized  hole 
in  a  disk  shaped  specimen  through  a  wedging  mechanism  was  developed,  but  stable  crack 
growth  could  not  be  achieved  in  the  stiffest  materials.  In  this  technique  the  load 
determination  from  the  wedging  crack  opening  mechanism  usually  suffered  from  load  train 
relaxations  at  Hertzian  contacts. 

Nevertheless,  an  apparently  significant  difference  in  toughness  in  materials  with  nominally 
similar  production  history,  was  found,  which  was  attributed  to  a  possible  microstructural 
parameter,  such  as  particle  size  distribution,  associated  with  the  highly  complex  three 
dimensional  microstructural  space  network.  In  the  approach  taken  here  to  isolate  the 
pertinent  elements  of  microstructure  design,  the  microstructure-property  connection  is 
studied  experimentally  with  new  or  improved  set  of  mechanical  property  testing 
procedures,  while  in  a  parallel  pair  of  computer  simulations  the  processes  that  control  the 
microstructure  [4]  and  the  mechanisms  that  control  the  resulting  fracture  response  [12]  are 
modeled.  Ideally,  this  approach  is  pursued  until  a  close  correspondence  is  reached  with  the 
actual  experimentally  probed  system  through  a  systematic  variation  of  the  available 
parameters  that  control  the  simulated  microstructures.  Once  a  proven  cause-effect 
relationship  is  reached  in  the  simulations  that  parallel  the  experimental  one,  it  becomes 
possible  in  principle,  to  use  the  computer  simulations  to  explore  other  microstructure- 
property  connections  and  determine  the  possible  benefits  of  pursuing  new  developments  in 
processing  to  achieve  such  microstructures.  This  present  communication  discusses  the 
results  obtained  for  the  connections  between  processing,  microstructure  and  mechanical 
properties  of  the  RBSN,  produced  by  Haggerty  and  co-workers  [6]. 


2.0  EXPERIMENTAL  DETAILS 


2.1  Production  of  Material 

Reaction  Bonded  Si3N4  samples  were  made  from  high  purity  Si  powder  synthesized  from 
CO2  laser  heated  SiH4  [13].  The  Si  powder  was  colloidally  dispersed  in  octanol,  pressed, 
dried  and  nitrided  at  1400C  in  UHP  N2.  Complete  procedures  for  production  have  been 
described  elsewhere  [14], 

2 . 2  Modulus  Measurement 

A  convenient  and  effective  way  of  measuring  the  average  elastic  properties  of  a 
heterogeneous  solid  of  a  small  microstructure  wave  length  is  by  acoustic  means  at  a  wave 
length  much  larger  than  the  microstructure  dimensions  [15,16].  Such  techniques  can  also 
be  used  to  detect  strength  limiting  flaws  of  dimensions  much  larger  than  the  wave  length  of 
the  acoustic  probe  [17-20]. 

In  order  to  probe  for  important  microstructural  differences  between  samples  the  Young’s 
Modulus  was  measured  by  an  ultrasonic  pulse-echo  technique  which  lends  itself  well  to  the 
measurement  of  elastic  properties  of  micro-porous  material  [21].  The  experimental 
approach  used  here  was  to  mount  a  thin  specimen  to  a  gold  coated  piezoelectric  single 
crystal  quartz  transducer  for  acoustic  transmission  as  depicted  in  Fig.  2.  The  optimal 
specimen  thickness  was  1mm,  which  is  thick  enough  in  relation  to  the  wavelength  of  the 

acoustic  pulse  (X  =  200|im),  and  thin  enough  so  that  minimal  attenuation  of  the  acoustic 
wave  occurred  and  a  plane  wave  condition  was  obeyed  [22].  The  volume  of  the  material 
sampled  was  approximately  S.lO'l^m^. 

The  acoustic  wave  was  generated  by  a  laser  pulse  of  1-1. 5mm  diameter  focused  on  the  gold 
film  which  creates  a  well  defined  pressure  pulse  that  is  substantially  planar  [22,23].  The 
pressure  pulse  forms  clear  signals  at  entry  into  the  sample  and  at  exit,  separated  in  time  and 
that  are  clearly  resolvable  by  the  piezo  electric  transducer.  In  this  way,  the  time  of  the 
transmission  of  the  stress  wave  into  the  specimen,  (and  out  of  the  transducer),  and  the  time 
of  return  of  the  echo  back  into  the  transducer  was  obtained.  The  difference  in  time  between 


these  two  events  is  the  time  of  travel  of  the  acoustic  wave  through  the  thickness,  reflection 
at  the  free  surface  and  return  to  the  specimen/transducer  interface.  The  uniaxial  strain  wave 
velocity  was  thus  obtained  and  is  related  to  the  Young’  Modulus  by: 

E  =  (l+T))(l-2a))pv2/(l-'0)  (1) 

where  p  is  the  material  density,  v  is  the  wave  velocity,  and  v  is  the  Poisson’s  ratio. 

2.3  Fracture  Toughness  Testing 

Fracture  strength  was  measured  using  the  biaxial  ball-on-ring  method  [24,25].  In  one 
method  the  fracture  toughness  was  related  to  the  strength  and  the  observed,  strength¬ 
controlling-flaw  size  [26],  by  curve  fitting  the  critical  flaw  size  and  site  stress  data  by 
examining  the  macroscopic  fracture  surfaces. 

In  a  new  and  more  accurate  experimental  technique  the  fracture  toughness  was  determined 
by  measuring  the  shape  of  the  open  crack  associated  with  a  Vickers  micro-hardness 
indentation.  The  motivation  for  developing  this  technique  was  to  obtain  fracture  toughness 
data  of  sufficient  resolution  to  make  a  more  precise  connection  with  microstructural 
characteristics.  Previous  attempts  to  measure  small  variations  in  fracture  toughness  in 
certain  RBSN  specimens  have  proven  unsuccessful  with  existing  indentation  crack  length 
methods  [27].  The  new  technique  could  be  shown  to  correspond  to  accepted  values  for 
reference  materials  for  which  toughness  values  are  known,  and  reveals  variations  in 
toughness  previously  too  small  to  measure.  This  technique  is  described  more  fully  in  a 
separate  communication  [28],  and  is  summarized  in  the  Appendix. 

2.4.  Microscopy 

Microscopy  of  the  RBSN  microstructure  involved  standard  grinding  and  lapping 
procedures  down  to  l|im  diamond,  which  produced  surfaces  with  many  filled-in  pores  that 

were  found  to  be  resistant  to  ultrasonic  cleaning.  However,  a  final  ion  plasma  etch  for  30 
minutes  at  6kV  removed  sufficient  material  to  reveal  the  microstructure  down  to  the 
submicron  pores.  Stereo-imaging  of  the  produced  surface  at  5000  magnification  revealed  a 


nearly  flat  surface  on  a  scale  larger  than  the  submicron  pores,  with  a  deviation  from  planar 
in  a  wave  like  form  with  a  wave  length  of  approximately  15  microns  and  an  amplitude  of 
0.2  microns.  This  long  wave  length  roughness  is  attributed  to  the  combination  of  lapping 
and  ion  milling  and  was  considered  acceptable  for  the  analyses  of  the  pore  features  which 
are  substantially  smaller  than  this  surface  wavelength. 

It  was  found  that  depositing  a  30nm  carbon  coating  allowed  for  good  high  voltage 
microscopy  and  resolution  of  sub-micron  features.  Due  to  the  low  atomic  mass  of  carbon, 
the  interaction  with  the  electron  beam  is  weak,  and  this  layer  is  mostly  penetrated  by  the 
high  energy  electron  beam  of  10-30kV  accelerating  voltage.  Objects  of  20nm  dimension  on 
such  ion  milled  surfaces  with  C  coating  have  been  resolved  with  a  high  resolution  field 
emission  SEM  (FESEM),  and  the  point-to-point  resolution  of  the  images  was  found  to  be 
approximately  5nm.  In  contrast,  for  uncoated  specimens  of  RBSN,  the  field  emission 
SEM  resolves  20nm  features  well  with  minimal  charging  at  2.0kV.  An  additional 
technique  used  for  observing  uncoated  specimen  surfaces  was  the  environmental  scanning 
electron  microscope  (ESEM).  This  microscope  operates  with  3  to  10  torr  of  water  vapor 
pressure,  which  prevents  excess  charge  buildup.  In  this  manner,  high  accelerating  voltages 
can  be  used  to  probe  non-conducting  specimens  without  standard  coatings.  Resolution 
obtainable  with  the  ESEM  is  approximately  20nm.  The  additional  feature  the  latter 
technique  provides  is  the  sampling  of  a  :.ub-surface  layer  of  material  into  which  the  high 
energy  electron  beam  penetrates  to  a  depth  that  is  a  function  of  the  material  properties.  The 
same  specimen  surface  imaged  in  the  ESEM,  without  and  with  a  20nm  Au  coating  is 
shown  in  Fig.  3.  The  uncoated  image  reveals  a  texturing  of  the  contrast  attributed  to 
density  fluctuations.  These  large  pore  channels  have  been  observed  previously  [29]. 

In'  order  to  objectively  determine  approximate  pore  sizes  and  shapes  in  highly  complex 
nano-scale  microstructures,  two  complementary  procedures  of  image  analyses  on 
metallurgically  prepared  surfaces  were  undertaken  in  parallel.  Stereo  images  with  SEM  of 
the  microstructure  in  the  first  case  were  hand  edited  to  differentiate  solid  from  pore.  In  the 
second  case  the  images  were  digitized  by  an  optical  scanner  and  pore  boundaries  were 
determined  by  a  contrast  threshold  criterion.  In  both  cases  the  images  were  statistically 
analyzed  with  a  numerical  image  analysis  system.  Pore  size  was  measured  by  pore  area 
and  the  dimensions  of  an  equal  area  ellipse.  The  shape  of  the  fitted  ellipse  was  determined 
in  such  a  way  that  maximum  overlap  occurred  with  the  pore  image.  The  pore  aspect  ratio 


was  characterized  by  the  ratio  of  the  major  and  minor  axes  of  the  equal  area  ellipse  and  the 
pore  surface  roughness  was  determined  by  the  measured  pore  circumference,  C, 

normalized  by  the  circumference  of  the  equivalent  area  circle,  C/2VTtA,  where  A  is  the  pore 
area.  The  circumference  was  measured  at  a  resolution  of  8.5nm. 

In  addition,  fracture  paths  were  studied  in  120  micron  thick  TEM  specimen  blanks  indented 
with  a  1kg  Vickers  indenter  before  impregnating  with  epoxy.  Microscopy  was  performed 
on  an  Akashi  EM002B  microscope.  More  detailed  description  of  the  microscopy  procedure 
is  given  elsewhere  [26]. 

2.5  X-ray  Diffraction  Experiments 

Ground  samples  of  9  micron  thickness  with  1  micron  finish  in  as-nitrided  states  were  used 
to  analyse  the  interior  and  exterior  regions  of  the  RBSN  pellets.  The  fraction  of  unreacted 
Si  the  alpha^eta  ratio,  grain  size,  lattice  parameter  and  lattice  distortion  were  determined  in 
this  manner  [26].  Lattice  distortion  can  originate  from  many  sources,  including  localized 
defects  and  microstrains;  they  should  be  considered  separate  from  changes  in  the  mean 
lattice  parameter  values,  which  result  from  long  range  macro-residual  stresses. 

3.0  EXPERIMENTAL  RESULTS 

3 . 1  Microstructure 

The  bulk  density  in  the  final  nitrided  samples  ranged  from  69.5%  to  83%  of  the  theoretical 
value.  The  percent  conversion  to  nitride  was  determined  by  x-ray  diffraction,  and  by 
weight  gain.  It  was  between  97-100%  [30],  and  was  uniform  throughout  the  specimens, 
within  2%.  The  grain  sizes  as  observed  by  TEM  and  x-ray  diffraction  were  between  45 
and  114  nm  respectively.  Lattice  distortion,  by  x-ray  line  broadening,  was  determined  to 
be  approximately  between  0.2  and  0.15%.  Processing  histories  for  specific  samples  can  be 
found  in  Table  2. 

Figure  4  shows  a  typical  microstructure  of  porous  RBSN  produced  by  Haggerty  and  co¬ 
workers  [6].  The  pores  appear  somewhat  equiaxed  with  relatively  smooth  surfaces 


indicating  that  during  the  formation  of  the  RBSN  there  is  a  finite  amount  of  time  for  some 
smoothing  of  the  surfaces  by  diffusion.  The  pore  size  distribution  in  the  typical 
microstructure  is  usually  narrow,  and  is  in  the  range  of  50-300nm,  as  shown  in  Fig.  5. 

Pore  aspect  ratios  were  found  to  range  between  1  and  4,  and  pore  circumferences 
normalized  with  that  of  the  equal  area  circle  range  between  1  and  2,  as  shown  in  Fig.  6. 
These  measures  reinforce  the  qualitative  observations  of  a  relatively  homogeneous, 
equiaxed  porosity  with  smooth  surfaces.  As  can  be  seen  in  Figs.  7a  and  7b  a  comparison 
of  pore  size  distributions  as  determined  by  the  major  axis  of  the  equal  area  ellipse,  between 
the  method  of  manually  edited  micrographs  and  that  obtained  by  a  digital  contrast  threshold 
method  revealed  no  significant  differences  in  average  dimension  and  distribution  of  pore 
size.  Manual  editing  produced  fewer  pores  with  diameter  less  than  50nm,  since  such  pores 

at  2x10^^  magnification  are  below  the  threshold  of  visual  detection.  All  measured  pore  size 
distributions  over  the  porosity  range  of  1 8  to  30%  were  alike  within  experimental  error, 
with  variations  arising  only  from  limited  sample  size,  as  seen  in  Figs.  5a  and  5b  for  the  two 
extremes  in  porosity.  In  a  few  instances  a  separate  population  of  pore  channels  was  found 
(see  Table  1).  These  pore  channels  shown  in  Fig.  8,  appear  as  high  aspect  ratio  features  of 
approximately  200nm  width  and  typically  more  than  1  micron  length,  and  form  an  apparent 
cell-like  network  with  a  periodic  length  scale  of  approximately  2-5microns.  In  addition,  the 
difference  in  area  fraction  of  micron  size  pores  observable  between  two  nominally 
identically  processed  materials:  i.e.,  one  with  the  elongated  pore  population  (specimen  G) 
and  the  other  with  quasi-homogeneously  distributed  equiaxed  pores  (specimen  H)  suggests 
that  the  volume  fraction  of  these  elongated  pore  channels  is  approximately  7%.  Line  profile 
height  measurements  of  the  sectioned  surface  revealed  that  5%  of  the  surface  sampled  was 
a  large  dimension  pore.  Profilometry  measurements  were  also  performed  by  a  scanning 

laser  confocal  microscope  which  has  a  height  resolution  of  0.25p,m  and  a  vertical  resolution 

of  l|i,m.  Qualitative  observation  suggests  that  at  least  some  of  the  elongated  pores  are 
interconnected,  providing  easy  fracture  paths. 

A  paragraph  discussing  determination  of  porosity  by  porosimetry  and  internal  pore  surface 
and  pore  connectivity  determination  by  gas  permeation  should  come  in  here. 


3.2  Modulus 


As  expected  the  modulus  was  inversely  dependent  on  the  average  porosity  as  is  shown  in 
Fig.  9.  In  the  case  where  the  microstructure  included  a  population  of  elongated  pores,  this 
tended  to  produce  a  lower  modulus,  due  to  the  interconnected  pore  channels.  The 
dependence  of  modulus  on  average  porosity  is  expected.  This  microstructural  morphology 
and  its  connection  with  the  modulus  over  the  range  of  porosity  mirrors  the  result  obtained 
from  the  quantitative  microstructure  analysis  described  above.  The  dependence  of  the 
measured  moduli  on  porosity  are  compared  in  Fig.  9  with  the  Eshelby  inclusion  model  of 
Chow  [31]  and  the  self  consistent  model  of  Budiansky  [32].  The  data  with  its 
characteristic  scatter  falls  somewhere  in  between  these  models.  The  degree  of  agreement  of 
the  data  with  the  models  reinforces  the  notion  of  a  quasi  random,  low  aspect  ratio  pore 
population  of  roughly  uniform  dimensions. 

Describing  the  effect  of  porosity  by  the  empirical  exponential  curve  fit  of  E  =  Eoexp(-bp) 
where  Eq  is  the  modulus  of  the  fully  dense  material  and  p  is  porosity,  requires  for  the  best 

fit  b  =  2.23  ±  0.17  as  compared  to  published  values  for  many  varieties  of  silicon  nitride  of 
2.4  ±  0.4  [33],  and  3.7  for  a  variety  of  reaction  sintered  silicon  nitrides  [34].  This 
relatively  low  value  of  b,  as  compared  to  the  range  of  values  of  1.7  to  6.6,  in  which  the 
large  majority  of  polycrystalline  materials  fall,  implies  that  for  the  isotropic  RBSN,  the 
elastic  properties  cannot  be  substantially  increased.* 

The  uncertainty  for  measured  values  of  E  with  this  technique  is  estimated  to  be  4%,  and  is 
associated  with  obtaining  a  uniform  thickness  specimen  and  measuring  this  dimension. 

3 . 3  Tensile  Strength 

Tensile  strengths  of  RBSN  with  porosity  ranging  from  17  to  27%  were  between  592  to 
308  Mpa,  respectively. 


*  A  lower  value  of  b  implies  less  dependence  on  the  averave  porosity 


3.4  Fracture  Toughness 


Fracture  toughness  determination  through  post-failure  fractographic  measurement  of 
strength  controlling  flaw  size  produced  fracture  toughness  values  in  the  range  of  1 .65  to 
2.7  MpaVm.  The  subjective  nature  of  determination  of  flaw  size  however,  limits  the 
accuracy  of  this  method. 

The  alternate  method  of  measurement  of  fracture  toughness,  by  determining  the  shapes  of 
open  cracks  under  load  [28]  produced  somewhat  smaller  values  of  fracture  toughness  when 
compared  with  values  obtained  by  the  method  based  on  flaw  size  measurements.  The 
values  ranged  from  0.93  to  1.97  MpaVm  and  scale  fairly  strongly  with  average  porosity,  as 
is  shown  in  Fig.  10.  The  few  cases  at  the  highest  levels  of  porosity  exhibited  toughness 
and  energy  release  rates  comparable  or  larger  than  values  of  material  with  substantially  less 
porosity  as  is  shown  in  Figs.  10  and  11.  These  highest  porosity  microstructures  exhibited 
evidence  of  increased  complexity  in  crack  path  in  comparison  with  the  largely  planar 
character  of  cracks  in  the  lower  porosity  materials.  The  fracture  path  micrograph  of  Fig.  12 
shows  some  evidence  of  crack,  bridging  behavior,  on  a  length  scale  of  the  order  of  the 
particle  diameter  in  the  green  body,  implying  some  connection  with  packing  defects. 

4.0  DISCUSSION 

The  microstructures  of  the  RBSN  presented  here  have  a  high  degree  of  uniformity  based  on 
microscopy  observations  and  the  relatively  small  scatter  in  the  elastic  modulus  data.  In 
addition,  the  similarity  of  pore  size  distributions  over  the  porosity  range  implies  a  degree  of 
qualitative  similarity  in  microstructural  features.  While  the  average  pore  diameter  is  related 
to  the  initial  Si  particle  diameter,  i.e.,  between  200  and  300nm,  the  pore  size  distribution 
based  on  stereological  measurement  performed  on  planar  sections,  shows  that  the 
microstructure  consists  of  many  very  small  pores  less  than  lOOnm.  These  small  pores, 
quite  likely,  correspond  to  bottleneck  configurations  between  larger  pores  and  appear  to 
have  little  consequence  on  the  overall  mechanical  properties,  as  evidenced  by  the  highly 
planar  character  of  the  typical  fracture  paths. 

Microstructural  homogeneity  was  found  frequently  except  in  the  few  cases  where  pores 
existed  on  a  scale  substantially  larger  than  the  average  intrinsic  pore  size.  The  processing 


steps  which  produced  these  defects  could  not  be  identified.  However,  the  length  scale  of 
the  pores  which  is  much  larger  than  tnc  initial  silicon  particle  size  suggests  that  the  source 
of  inhomogeneity  is  associated  with  particle  packing  defects  or  stress  relief  cracks,  and  may 
provide  a  low  resistance  crack  propagation  path.  These  materials  indeed  exhibited  lower 
strength  and  toughness. 

The  grain  size,  as  determined  in  two  representative  samples,  was  found  to  be 
approximately  one  half  of  the  average  pore  diameter,  thus  grain  boundary  surface  area  is 
comparable  with  pore  surface  area  at  these  relatively  high  porosities.  This  result  is 
consistent  with  the  fact  that  the  fracture  paths  were  observed  to  be  predominantly 
intergranular  from  observations  of  indentation  cracks  viewed  by  TEM.  This  could  be  due 
to  either  weaker  grain  boundaries  or  due  to  grain  level  internal  stresses  of  misfit.  The  fact 
that  the  latter  may  be  important  and  that  they  may  have  a  wavelength  larger  than  the  grain 
size  might  be  implied  by  the  apparent  stress  relief  cracks  observable  in  some 
microstructures. 

The  monotonic  decrease  of  elastic  modulus  with  average  porosity  is  not  fully  paralleled  by 
the  decrease  of  fracture  toughness  and  energy  release  rate.  The  highly  planar  character  of 
the  fracture  paths,  especially  at  low  porosity,  highlights  the  small  consequence  of  the  pore 
configurations,  and  the  relative  importance  of  the  average  porosity  level.  At  the  highest 
porosity  levels  (p  =  30%),  the  observed  increased  crack  complexity  and  bridging  differs 
markedly  from  a  planar  section  probe  of  the  microstructure,  and  is  paralleled  by  an 
increased  fracture  toughness. 

The  simulation  of  fracture  paths  in  similar  porous  microstructures  [12]  illustrates  some  of 
the  results  already  outlined.  Highly  uniform  microstructures  have  little  potential  for 
harnessing  toughening  mechanisms  such  as  crack  tip  shieldings  by  redundant  cracking, 
crack  bridging  or  crack  deflection.  However,  in  such  simulations  it  was  found  that 
increasing  local  material  variability  in  terms  of  strength  and  stiffness  variation  in  the  porous 
material  on  the  scale  of  the  typical  distances  of  coherent  crack  advance  such  as  the  grain 
size,  which  was  found  to  increase  the  crack  complexity  and  associated  toughness.  In 
addition,  increased  microstructural  variability  in  terms  of  strength  and  stiffness  was  found 
to  arise  from  increasing  amounts  of  porosity  on  length  scales  equivalent  to  a  small  multiple 
of  the  pore  dimensions.  The  2-D  porosity  level  for  optimal  strength  variability  occurs 


within  the  range  of  10-30%.  Less  porosity  corresponds  to  low  material  variability  and  high 
porosity  implies  high  variability  and  overall  material  weakening.  However,  in  view  of  the 
fact  that  material  strength  and  stiffness  correlate  strongly**  and  was  found  to  significantly 
reduce  crack  complexity,  it  is  likely  that  increased  crack  complexity  may  be  found  at 
porosity  levels  higher  than  those  observed  here. 

New  directions  for  further  improvement  of  toughness  include  investigating  the  specific 
details  of  the  increase  in  crack  complexity  at  high  porosity  levels.  Clearly,  at  levels  of 
quasi-random  porosity  less  than  the  3-D  elastic  percolation  limit,  small  interaction  between 
pores  will  generate  fracture  behavior  appropriate  for  dilute  porosity  such  as  planar  fracture 
and  average  sampling  of  material.  At  the  highest  levels  of  porosity  studied  here,  when 
porosity  can  no  longer  be  randomly  distributed  the  crack  complexity  was  found  to  be  only 
marginally  greater  than  that  in  material  with  planar  fracture  paths.  Therefore,  a  further 
increase  in  porosity  with  a  qualitatively  similar  microstructure  morphology,  may  reveal 
significantly  higher  amounts  of  crack  complexity,  such  as,  roughness  and  bridging.  The 
increase  in  crack  tip  shielding  effects  at  higher  porosities,  however,  must  out-weigh  the 
reduction  in  available  solid  fracture  surface  area  and  average  elastic  modulus.  These  trends 
illustrate  the  competition  between  the  increase  in  toughness  due  to  higher  variability  in 
strength  and  stiffness  at  higher  porosity  and  the  reduction  in  toughness  due  to  the  reduction 
of  the  elastic  and  strength  properties  with  increasing  porosity. 


'*  A  control  microstructure  in  the  simulation  deemed  to  have  characteristics  of  the  real  RBSN  exhibited 
strength  and  stiffness  variation  of  material  elements  on  the  scale  of  a  multiple  of  the  pore  size.  These 
variations  were  correlated  by  the  two  variable  correlation  parameter,  x  =  0.85.  Although  a  high  correlation 
of  strength  and  stiffness  through  porosity,  for  crack  tip  material  elements  was  generally  found,  the 
dependece  was  weakened  at  higher  porosities.  The  implication  is  that  increased  material  variability  at  higher 
porosity  produces  a  de-coupling  of  local  strength  and  stiffness,  and  results  in  comparatively  higher  crack 
complexity  and  toughness. 


Other  possible  avenues  of  potential  toughness  improvement,  as  suggested  by  the  above 
mentioned  simulation  results,  involve  exploring  higher  variability  microstructures.  Material 
variability  contributions  were  found  to  .  nse  in  simulated  fracture  with  the  introduction  of 
grain  boundary  strength  less  than  the  bulk  matrix  strength,  suggesting  grain  boundary 
strength  tailoring.  Other  potential  sources  of  increased  material  variability  include  thermal 
mismatch  residual  stresses  and  elastic  anisotropies  which  are  considered  as  effective 
heterogeneities  in  the  nucleation  of  microcracks  stimulated  by  the  main  crack.  Other  typical 
toughening  mechanisms  such  as  crack  wake  bridging  suggest  higher  grain  aspect  ratio  with 
random  principal  directions  and  the  incorporation  of  whiskers. 

5.0  CONCLUSIONS 

The  high  degree  of  homogeneity  of  microstructures  of  RBSN  on  scales  larger  than  the 
average  sub-micron  pores  and  the  generally  planar  character  of  the  observed  fracture  paths 
implies  that  the  likelihood  of  obtaining  significantly  improved  fracture  toughness  is  not 
high  in  material  with  such  quasi-homogeneous  microstructures.  At  the  highest  levels  of 
porosity  observed  in  this  study  (30%)  increased  crack  complexity  and  toughness  was  in 
evidence.  Higher  levels  of  porosity  have  the  potential  for  increased  levels  of  fracture 
energy,  as  has  been  found  by  Rice  et  al  [7],  and  was  noted  in  fracture  path  simulations 
[12].  Higher  porosity  and  other  factors  such  as  grain  boundary  strength  tailoring,  was 
found  to  provide  increased  variability  and  increased  crack  complexity  and  toughness. 
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Appendix  Fracture  Toughness  Determination  through  Measurement  of 
Open  Crack  Spaces. 

The  stress  state  at  the  crack  tip  of  Vickers  indentation  induced  cracks  is  obtainable  from  the 
gradient  of  separation  between  the  crack  faces  as  a  function  of  distance  from  the  crack  tip 
[28].  Thus,  the  accuracy  of  this  method  relies  on  the  fundamental  connection  between 
observable  open  crack  shape  on  the  verge  of  propagation  and  the  critical  crack  tip  stress 
intensity.  Through  relating  computed  crack  shapes  obtained  by  finite  element  solutions  for 
boundary  loading  conditions  appropriate  to  the  indentation  process,  to  open  crack  shapes 
determined  by  SEM,  it  was  possible  to  obtain  high  quality  data  for  the  fracture  toughness. 

Figure  13  shows  the  open  crack  displacements  Av  for  a  specific  case  of  a  Vickers 

indentation  crack  along  the  length  of  the  crack,  fitted  against  a  matching  finite  element 
solution  for  a  control  experiment  on  soda  lime  glass.  The  fracture  toughness  determined  by 

the  fit  was  Kjc  =  0.55MPaVm. 

The  implementation  of  the  measurement  technique  involves  the  preparation  of  the  surface 
for  microscopy,  the  indentation  of  the  surface  with  the  Vickers  microhardness  technique 
[35]  and  generation  of  high  resolution  micrographs  of  the  shapes  of  the  produced  cracks. 

Prior  to  indentation,  a  standard  gold  coating  used  regularly  in  SEM,  is  applied  for  the 
purpose  of  obtaining  optimal  resolution.  Indentations  are  made  with  a  standard  Vickers 
microhardness  apparatus  with  loads  ranging  from  2  to  ION.  The  choice  of  load  is 
important  to  the  generation  of  radial  cracks  which  are  large  enough  for  observation. 
However,  care  must  be  taken  not  to  exceed  the  load  for  which  a  profusion  of  other  cracks 
are  produced.  The  load  is  held  for  20  seconds,  and  removed.  Microscopy  was  performed 
within  several  minutes  of  indentation,  and  with  sufficient  resolution  to  clearly  delineate  the 
crack  edges. 

Open  crack  shapes  are  obtained  directly  from  the  micrographs  with  standard  micrometer 
measurements.  Care  was  taken  to  estimate  the  uncertainty  of  open  crack  flank  separation 
measurement  by  providing  upper  and  lower  bounds,  as  observed  from  the  micrographs. 
The  upper  and  lower  bounds  were  determined  by  the  uniformity  of  the  contrast  gradients. 
The  technique  pre-supposes  that  the  indentation  crack  is  kept  open  by  fracture  debris 
pressed  into  the  crack  by  the  indent  and  that  crack  closing  is  negligible. 


Nevertheless,  due  to  this  assumption  the  reported  fracture  toughness  values  are  considered 
to  be  under  estimates. 


V 


Fig. 2.  The  pulse-echo  method  of  modulus  measurement  by  means  of  a  piezo  electric 
transducer. 


Fig. 3.  McUciuil  C  obscived  in  the  SEM  v\'ith  an  iincoated  sui'lace  (aj;  and  with  a  surlace 
coated  with  a  20nm  Au  film  (b). 


Fig. 4.  Typical  micro.staicture  of  RBSN  showing  quasi-uniformly  distributed  submicron 
pores. 
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Fig. 5.  Two  typical  examples  of  pore  size  distribution  in  materials  on  the  extremes  of  the 
porosity  range:  (a)  material  M  with  18%  porosity:  (b)  material  A  with  30.5% 
porosity. 
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Fig. 6.  Typical  distributions  of  pore  aspect  ratio  (a);  and  pore  circumference  normalized 
with  the  circumference  of  a  circle  with  equal  area.  (b). 
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Fig. 7.  Comparison  of  methods  for  determination  of  pore  characteristics:  (a)  porosity 
characterized  from  manually  edited  micrographs;  (b)  porosity  characterized  by 
automated  image  analyses. 
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Fig. 9.  Effect  of  porosity  on  modulus  is  compared  with  two  theoretical  models:  the 

Eshelby  inclusion  model  of  Chow  [31]  and  the  self  consistent  model  of  Budiansky 
[32]. 
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Fig.lO.Dependence  of  fracture  toughri'-  ;s  of  RBSN  on  porosity.  The  range  of  values  for 
each  material  represents  the  local  toughness  variations  and  do  not  represent  error 
bars  of  macro  tests. 
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Fig.  13.  A  specific  test  case  fit  for  the  determination  of  fracture  toughness  from  open  crack 


shape  measurements  in  soda  lime  glass  where  a  fracture  toughness  of  Kjc  =  0.55 


Mpa  Vm  was  determined  from  the  fit  of  experimental  results  to  the  crack  shape 


obtained  by  finite  element  computation. 
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A  NEW  METHOD  OF  FRACTURE  TOUGHNESS  DETERMINATION 
IN  BRITTLE  CERAMICS  BY  OPEN  CRACK  SHAPE  ANALYSIS 


F.  Haubensak*  and  A.S.  Argon^ 

Massachusetts  Institute  of  Technology,  Cambridge,  MA  02139 


ABSTRACT 

Improvement  of  the  fracture  toughness  of  high  quality  ceramics  remains  as  one  of  the  most 
important  goals  in  materials  development.  An  associated  problem  is  the  accurate 
measurement  of  fracture  toughness  in  such  brittle  or  semi-brittle  ceramics,  particularly  in 
small  samples  encountered  in  material  development.  Previously  used  methods  relying  on 
measurement  of  the  size  of  fracture  mirrors,  the  indentation  load  and  crack  length  in 
Vickers  hardness  induced  cracking,  and  a  variant  of  similar  techniques  have  all  been  less 
than  satisfactory  in  discriminating  quantitative  differences  among  materials.  A  hitherto 
unused  technique  of  inferring  the  fracture  toughness  in  samples  from  measurements  of 
open  crack  flank  displacements,  vv'hich  we  have  developed,  avoids  most  of  the  theoretical 
and  experimental  difficulties  of  other  methods.  While  it  is  somewhat  intensive  in  terms  of 
evaluation  and  requires  high  resolution  of  open  cracks,  the  technique  is  fundamentally  the 
soundest  of  all  techniques  and  is  capable  of  furnishing  discriminating  results.  We  present 
results  of  its  application  to  the  measurement  of  some  model  materials  such  as  soda  lime 
glass,  single  crystal  silicon,  alumina,  and  a  reaction  bonded  silicon  nitride  whose  porosity 
would  ordinarily  present  difficulties  with  other  methods. 


*  Current  address:  Department  of  Materials  Science  and  Engineering,  Stanford  University,  Stanford,  CA 
94305-2205 

^  Author  for  correspondence. 


1.0  INTRODUCTION 


In  the  development  of  materials  with  superior  fracture  toughness  the  necessity  of  high 
resolution  in  the  toughness  measurements  on  small,  laboratory  scale  specimens  is 
paramount.  In  this  pursuit,  a  number  of  methods  have  been  applied  in  the  development  of 
high  quality  reaction  bonded  porous  silicon  nitride  (RBSN)  produced  by  Haggerty  et.  al. 

[1] .  Nominally  defect-free  disk  specimens  initially  were  fractured  in  a  ball-on-ring  test 
configuration,  and  the  fracture  mirrors  were  measured  and  related  to  fracture  toughness 
through  the  bend  strength.  These  results  showed  broad  trends,  albeit  with  large  uncertainty 

[2] .  Standard  indentation-toughness  techniques,  well  described  in  the  literature  [3,4],  were 
also  used,  but  proved  to  lack  sensitivity.  An  experimental  technique  to  pry  open  cracks 
emanating  from  circular  holes  in  plate  specimens  produced  excellent  results  in  model 
materials,  but  was  limited  in  use  with  the  high  stiffness  silicon  nitride  material  of  interest. 
Large  corrections  in  the  induced  displacements,  due  to  compression  of  the  wedging 
mechanism,  while  calculable,  became  larger  than  the  actual  displacements,  with  the  result 
that  inaccuracies  became  too  large  to  be  acceptable. 

Finally,  it  is  known  that  the  crack  tip  stress  field  can  be  deduced  from  the  boundary 
conditions  and  the  crack  flank  opening  displacements,  if  the  elastic  properties  are  known. 
Therefore,  this  methodology  was  developed  further  and  for  the  special  boundary  value  case 
of  the  internally  pressurized  half  penny  shaped  surface  crack  associated  with  the  radial 
cracks  of  Vickers  indentations.  The  development  of  this  technique  and  its  use  in  measuring 
Kjc  in  RBSN  is  the  subject  of  the  present  communication. 

2.0  METHODOLOGY 

2. 1  Theoretical  preliminaries. 

In  a  Mode  1  crack  tip  field  for  a  crack  under  impending  motion,  Kj  =  Kjc  can  be 
determined  from: 


(1) 


kicJCQDI^e 

F(r/a)Va 

where  F(r/a)  is  a  non-dimensional  characteristic  crack  shape  factor  dependent  on  the 
geometry  and  the  mode  of  loading,  a  is  the  crack  length,  and  E  is  the  Young’s  modulus. 
The  function  can  be  determined  by  FEM  methods  for  a  crack  of  a  given  geometry  and  mode 

of  loading  under  a  unit  K^,  Va  and  E.  An  example  is  presented  in  Fig.  la  for  the  case  of  a 

Vickers  indentation  of  interest  in  an  elastic,  isotropic  material.  It  was  found  that  the  shape 
of  the  open  cracks  qualitatively  fit  best  with  the  solution  for  an  internally  loaded  crack  with 
point  loads  at  the  positions  corresponding  to  the  comers  of  the  indentation.  This  can  be 
seen  in  Fig.  lb,  where  various  open  crack  COD  solutions  are  compared  with  measurements 
from  indentations  produced  in  a  glass. 

A  fundamental  assumption  here  is  that  indentation  of  material  misfit  or  comminution  chips 
under  the  indentor  keep  the  crack  open  at  the  shape  for  its  impending  propagation.  There 
are  reasons  to  believe  this  to  be  the  case  in  many  circumstances  [5-7],  and  various 
experimental  observations  [8],  and  analytical  [9]  modeling  support  the  notion  that  the 
residual  stress  state  upon  the  removal  of  the  indentation  load  is  sufficient  to  extend  the 
radial  cracks  during  the  unloading  process  to  conditions  of  crack  arrest. 

Nevertheless,  because  of  such  possible  relaxations  from  a  critically  stressed  crack  state  the 
determined  fracture  toughness  values  must  be  regarded  as  low  toughness  estimates. 


2.2  Method  of  Observation 


Determination  of  the  crack  length  and  opening  shape  was  done  primarily  by  standard 
scanning  electron  microscopy  (SEM)  techniques,  in  order  to  attain  high  enough  resolution, 
typically  up  to  2x10^  magnification.  However,  other  methods  such  as  atomic  force 
microscopy  and  sometimes  even  light  optical  microscopy,  have  the  potential  to  generate  a 
similar  quality  of  data. 

Elastic  properties  were  obtained  either  from  the  literature  for  the  model  materials,  or  were 
measured  with  a  pulse-echo  technique  described  elsewhere  [10]. 

2.3  Measurement  of  Open  Crack  Shapes  from  Micrographs 

While  other  methods  exist  for  determining  crack  flank  opening  displacements,  a  standard 
SEM  was  utilized  for  the  measurements  obtained  here.  The  relative  opening  displacements 
were  taken  manually  from  the  micrographs.  Since  the  edges  of  the  crack  faces  are  not 
perfectly  defined,  upper  and  lower  limits  were  estimated  in  the  following  manner.  An 
upper  limit  was  determined  based  on  the  notion  that  beyond  a  certain  dimension,  the 
contrast  of  the  image  does  not  change  and  represents  the  contrast  of  the  surrounding  solid. 
In  the  same  way  the  estimate  of  the  lower  limit  was  determined  by  considering  the  contrast 
inside  the  crack. 

The  resolution  of  the  measurements  is  fundamentally  limited  by  the  resolution  of  the 
micrographs.  However,  the  capability  of  measuring  multiple  crack  opening  displacements 
permits  this  to  be  improved.  Obtaining  a  well  documented  shape  of  a  single  crack  reduces 
the  errors  in  comparison  with  a  single  measurement.  Typical  uncertainties  associated  with 
the  measurement  of  the  fracture  toughness  obtained  from  one  crack,  as  defined  by  the 
standard  error  parameter  was  found  to  be  3%. 


2.4  Results  of  Calibration  Experiments  in  Standard  Materials. 


In  order  to  determine  the  overall  accuracy  of  the  method,  the  measurement  procedure  was 
applied  to  three  model  materials:  silica  glass,  an  alumina  single  crystal  and  a  silicon  single 
crystal.  The  crack  opening  displacement  data  are  shown  in  Figs.  2,3,  and  4.  The 
corresponding  fracture  toughness  values  are  plotted  in  Fig.  5  for  the  data  from  the  glass 
with  the  above  described  procedure,  and  shows  an  independence  of  the  toughness  with 
position  along  the  crack.  This  implies  a  good  match  between  the  calculated  and  measured 
shape  of  the  open  crack  dimensions. 

Furthermore,  the  numerical  values  of  fracture  toughness  correspond  well  with  the  accepted 
values  for  the  three  model  materials,  as  shown  in  Table  1. 

3.0  RESULTS 

Figures  6  and  7  show  micrographs  of  the  porous  RBSN  with  indentation  induced  cracks, 
demonstrating  the  open  crack  configurations.  Figure  7  shows  an  example  of  the  frequently 
observed  bridging  phenomenon  of  a  ligament  across  the  crack  faces. 

These  indentations  in  the  RBSN  were  produced  by  Vickers  hardness  indenters  with  400g 
loads,  in  materials  which  ranged  in  porosity  from  10  to  30%.  In  each  case,  cracks  were 
consistently  produced  with  well  formed  shapes,  i.e.,  four  straight  radial  cracks  emanating 
at  right  angles  from  neighboring  cracks. 

The  calculated  fracture  toughness  values  as  a  function  of  porosity  are  shown  in  Fig.  8.  As 
has  been  found  in  other  porous  material  systems,  fracture  toughness  values  depend  on  the 
average  porosity  in  an  exponential  fashion  [11].  If  the  fracture  toughness  is  a  function  of 
the  porosity: 


(2) 


K  =  K  exp(-bp) 
o 

where  p  is  the  average  porosity,  then  b  is  found  to  be  2.3,  and  compares  well  with  data 
found  in  the  literature  where  b=2.4  a  range  of  silicon  nitrides.  Energy  release  rates  show 
this  trend  more  clearly  in  Fig.  9.  However,  at  the  high  porosity  extreme,  the  toughness 
values  tended  to  deviate  considerably  from  the  overall  trend.  This  deviation  correlates  with 
the  observation  of  a  discernible  increase  in  crack  tortuosity  in  the  highest  porosity 
materials,  as  seen  in  Fig.  7. 

The  two  exceptions  to  the  overall  trend  are  due  to  distinctly  different  processing  steps.  In 
material  (a)  a  binder  was  added  that  increased  the  final  density,  and  was  postulated  to  create 
a  final  microstructure  with  more  integrity.  Material  (b)  exhibited  large  scale  channel  shaped 
pores  which  provided  preferential  fracture  paths.  Nominally,  this  material  was  produced 
with  a  slightly  higher  concentration  of  fine  agglomerates  which  reduced  the  overall  material 
homogeneity. 

4.0  DISCUSSION  AND  CONCLUSIONS 

In  the  present  experimental  study  we  have  utilized  a  theoretically  well  known  tool  for 
determining  the  critical  stress  intensity  for  crack  growth  from  the  measurement  of  open 
crack  shapes.  The  method  has  been  shown  to  possess  good  accuracy  with  some  chosen 
model  materials  with  well  defined  toughness  values  and  has  higher  resolution  than  some 
previous  methods  involving  measurement  of  indentation  crack  lengths  associated  with 
indentation  load  measurements.  In  the  latter  determinations  of  fracture  toughness  based  on 
indentation  loads,  large  uncertainties  often  exist  resulting  from  frictional  resistances  which 
have  no  role  in  the  propagation  of  the  indentation  induced  cracks.  In  comparison,  the  crack 
shape  measurement  employed  by  us  avoids  such  pit  falls  by  relying  entirely  on  the  direct 


measurement  of  relative  displacement  gradients  along  the  crack  which  are  linked  to  the 
critical  crack  tip  stress  intensity.  The  clear  advantage  of  the  method  which  has  uncertainties 
less  than  3%,  are  to  some  degree  offset  by  the  increased  effort  required  in  preparing  high 
resolution  micrographs. 

The  method  of  measuring  through  measurement  of  crack  shape  was  utilized  here  to 

determine  influences  of  specific  processing  steps  on  the  fracture  toughness  of  a  series  of 
porous  silicon  nitrides  based  primarily  on  their  relative  porosities.  While  other  methods  of 
measuring  the  fracture  toughness  on  small  size  RBSN  specimens  indicated  negligible 
differences,  small,  but  definite,  fracture  toughness  differences  could  be  reliably  resolved 
with  the  new  technique.  The  measurements  established  a  definite  trend  of  decreasing 
toughness  with  increasing  average  porosity  which  has  been  commonly  found  in  other 
studies  as  well  [11].  However,  at  the  highest  levels  of  porosity  where  high  crack  tortuosity 
was  found,  fracture  toughness  values  deviated  significantly  upward  from  this  decreasing 
trend.  This  effect  appeared  to  be  due  to  increased  occurrence  of  crack  deflection  and  of 
crack  face  bridging  ligaments  in  the  high  porosity  materials. 
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Table  1.  Values  of  for  Model  Materials 


Material 

Measured  Fracture 

Previously  Reported 

Toughness 

Fracture  Toughness 

Soda  Lime  Glass 

0.81  ±  0.025 

0.40  -  0.85 

Si  -  single  crystal  [110] 

0.84  ±  0.03 

0.87 

a- Alumina  single  crystal 

2.86  ±0.10 

2.93* 

[110] 

value  for  hep  alumina 


6.00 


Figure  la.  The  non-dimensional  open  crack  displacement  function  for  an  isotropic  elastic 
2-D  material.  In  this  calculation,  the  application  of  the  point  load  dipole  occurs  at 
r/a  =  0.65. 


000 


Figure  lb.  Comparison  of  various  open  crack  shape  solutions  for  internally  loaded  surface 
cracks.  Solution  (1)  is  a  center  loaded  point  force,  (2)  is  for  uniform  internal  pressure  over 
the  portion  of  the  crack  representing  the  indentation  impression  (0.65a<r<a),  (3)  is  for 
symmetric  point  loading  about  the  impression  at  r/a  =  0.65,  similarly  (4)  is  for  point 
loading  at  r/a  =  0.40,  and  (5)  is  for  a  centrally  loaded  3-D  penny  crack. 
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Figure  2.  Open  crack  displacement  data  from  indentation  induced  cracks  in  silica  glass. 
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Figure  3.  Open  crack  displacement  data  from  indentation  induced  cracks  in  single  crystal 
alpha  alumina. 


o 

o 

o 

o 

o 

o 

in 

o 

in 

o 

in 

CM 

CM 

T“ 

T“ 

[LUU]  ‘^U0LU0OB|dSjQ  6U!U0dO  >10BJO 


Figure  4.  Open  crack  displacement  data  from  indentation  induced  cracks  in  single  crystal 
silicon. 


Figure  6.  Micrograph  of  mdeniation  uiduced  radial  crack  iti  RliSN. 


exhibiting  high  crack  dctlection  and  cvidcr  rc  It'--  bridging  ligaments. 


Figure  8.  Dependence  of  fracture  toughness  on  porosity  in  RBSN. 
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Figure  9.  Dependence  of  energy  release  rate  on  porosity  in  RBSN, 
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APPENDIX  D.  “A  Simulation  of  the  Reaction  Bonding  of  Si3N4 
to  Generate  Porous  Microstructures” 

F.G.  Haubensak,  V.V.  Bulatov,  and  A.S.  Argon, 
submitted  to  J.  Computer  Aided  Materials  Design. 


A  SIMULATION  OF  THE  REACTION  BONDING  OF  Si^N^TO  GENERATE 

POROUS  MICROSTRUCTURES 

F.  G.  Haubensak,*  V.V.  Bulatov,  and  A.S.  Argon** 

Massachusetts  Institute  of  Technology,  Cambridge,  MA,  02139 

ABSTRACT 

A  stochastic  numerical  technique  is  presented  that  has  been  used  to  study  the  formation 
of  reaction  bonded  porous  silicon  nitride  (RBSN),  Like  many  reaction  bonded  materials, 
RBSN  is  a  structural  ceramic  which  has  a  naturally  occurring  porosity  of  approximately  25% 
resulting  from  the  unique  kinetics  of  formation.  Since  such  porosity  is  usually  unavoidable, 
it  is  of  interest  to  discover  its  relevant  measures,  such  as  the  size  distribution  and  degree 
of  homogeneity  of  pores,  and  the  factors  that  govern  the  formation  of  these,  in  order  to 
understand  and  optimize  the  crack  propagation  resistance  of  such  materials.  The  probing  of 
the  fracture  toughness  of  the  generated  microstructures  by  computational  means  is  described 
in  the  companion  communication. 

The  primary  purpose  of  the  simulated  reaction  was  the  generation  of  realistic  microstruc¬ 
tures  incorporating  the  major  features  of  reaction  bonding.  In  this  manner  it  is  possible  to 
generate  objective,  systematic  variations  of  similar  microstructures.  With  such  microstruc¬ 
tures,  crack  nucleation  and  propagation  mechanisms  can  be  observed  and  processing  param¬ 
eters  which  control  fracture  toughness  can  be  identified. 

In  the  present  study,  we  outline  the  details  of  the  reaction  kinetics  model,  demonstrate 
its  ability  to  generate  realistic  microstructures,  and  report  some  additional  observations  on 
some  topological  phenomena  resulting  from  the  reaction,  such  as  those  which  control  the 
formation  of  large  pores  in  the  final  structure.  In  yet  a  separate  communication  of  a  related 
experimental  study  we  present  the  microstructures  of  the  actual  RBSN.  We  will  use  the 
findings  of  the  latter  here  for  making  some  comparisons. 


*  present  address,  Department  of  Materials  Science,  Stanford  University,  Stanford  CA 

**  author  for  correspondence 
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1.0  INTRODUCTION 


The  cause-effect  relationship  between  microstructure  and  mechanical  properties  of  ma¬ 
terials,  particularly  their  fracture  toughness  has  been  of  interest  for  some  decades.  In  most 
instances,  a  material  with  a  given  microstructure  is  studied  by  probing  it  with  a  number  of 
more  or  less  conventional  tests.  The  results  are  then  related  and  justified  through  a  mech¬ 
anistic  model  based  on  micromechanical  analysis.  In  a  few  instances  the  findings  of  the 
microstructure-property  connection  is  used  to  create  different  microstructures  which  should 
have  better  properties.  In  a  complementary  approach,  that  is  usually  only  taken  on  paper, 
more  ideal  microstructures  are  considered  which  through  the  existing  models  should  lead  to 
superior  performance.  The  usual  constraint  on  this  latter  approach  is  the  inadequacy  of  the 
available  degrees  of  freedom  in  processing  to  achieve  the  sought  after  ideal  microstructures. 
In  a  hybrid  approach  which  we  have  taken,  the  microstructure-property  connection  is  stud¬ 
ied  experimentally  while  in  a  parallel  computer  simulation  the  processes  that  control  the 
microstructure  and  the  mechanisms  that  control  the  property  are  modeled.  This  is  pursued 
until  a  close  correspondence  is  reached  with  the  actual  experimentally  probed  system  through 
a  systematic  variation  of  the  parameters  that  control  the  simulated  microstructures.  Once 
a  proven  cause-effect  relationship  is  reached  in  the  simulations  that  parallels  the  experimen¬ 
tal  one,  it  becomes  possible  to  use  the  computer  simulations  to  explore  other  conceivable 
microstructure-property  connections  and  determine  the  possible  benefits  of  pursuing  new 
developments  in  processing  to  achieve  such  microstructures. 

In  the  present  communication  we  develop  a  numerical  model  of  the  gas  phase  reaction 
that  results  in  porous  RBSN  for  the  principal  purpose  of  generating  non-arbitrary  microstruc¬ 
tures  that  are  topologically  similar  to  those  of  the  real  RBSN,  that  have  been  reported  by 
Haggerty  and  co-workers  [1-3].  In  a  companion  communication  we  have  explored  the  fracture 
resistance  of  these  computer  generated  microstructures  by  means  of  a  separate  simulation  of 
tensile  fracture  and  of  crack  propagation  through  these  porous  microstructures  to  develop 
the  microstructure-property  connection  [4].  This  pair  of  simulations  has  given  a  detailed 
and  valuable  understanding  of  the  gaseous  reaction  that  forms  the  porous  RBSN,  and  the 
mechanical  worth  of  the  specific  porous  microstructures. 

Figure  1  shows  a  typical  microstructure  of  porous  RBSN  as  produced  in  the  Ceramics 
Processing  Research  Laboratory  at  MIT.  The  pores  are  equiaxed  with  relatively  smooth 
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surfaces  indicating  that  during  the  formation  of  the  RBSN  there  is  a  finite  amount  of  time 
for  some  smoothing  of  the  surfaces  by  diffusion.  The  phase  content  and  grain  size,  as 
determined  by  TEM  and  x-ray  techniques  were  found  to  be  predominantly  of  the  a  phase 
at  75-80%  which  in  some  cases  tended  to  form  whisker-like  shapes,  while  the  remainder 
was  /3  Si^Ni  in  typically  larger  and  rod-like  morphology.  The  TEM  observations  indicated, 
however,  more  equiaxed  grains  of  a  size  roughly  equal  to  the  size  of  the  small  pores  in  the 
final  reacted  microstructure.  Moreover,  a  small  residual  concentration  of  unreacted  silicon 
was  also  detectable  which,  however,  is  difficult  to  capture  on  micrographs.  The  experimental 
observations  of  fracture  paths  established  that  fractures  follow  predominantly  along  grain 
boundaries  [3]. 

In  the  present  pursuit,  the  process  of  interest  has  been  the  reaction  bonding  of  silicon 
nitride  (RBSN).  In  a  closely  related  experimental  program,  Haggerty  and  co-workers  [1,2] 
have  produced  a  high  quality  porous  RBSN  from  sub-micron  size  spherical  Si  powder.  The 
microstructure  and  mechanical  properties  of  this  RBSN,  including  its  fracture  toughness, 
have  been  studied  in  considerable  detail  [3,5].  An  intriguing  aspect  of  these  experiments  was 
the  measured  fracture  toughness  which  ranged  from  1.0  to  2.7  MPa\/m  and  appeared  to  be 
related  to  different  and  reproducible  features  of  the  microstructure  [3].  Since  the  experimen¬ 
tally  observed  fracture  paths  showed  some  evidence  of  crack  deflection,  crack  bridging  and 
possibly  even  crack  tip  shielding,  resulting  from  some  surrounding  redundant  microcrack¬ 
ing,  the  possibility  of  improving  the  fracture  toughness  further  by  deliberate  microstructure 
manipulation  became  of  interest.  Such  a  pursuit  becomes  possible  in  the  pair  of  computer 
simulations  of  microstructure  and  crack  propagation,  presented  here  and  in  the  companion 
study  [4]. 
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2.  METHODOLOGY 


2.1  Reaction  Chemistry 

In  order  to  generate  realistic  porous  microstructures  of  RBSN,  it  is  important  to  include 
all  the  major  features  of  consequence  of  the  reaction  bonding  process  in  the  simulation.  The 
reaction  begins  with  the  formation  of  discrete  Si3N4  nuclei  which  become  established  over 
an  initial  incubation  period  [1].  These  nuclei  grow  at  the  expense  of  the  evaporating  silicon 
and  form  individual  grains  in  the  final  microstructure.  Under  typical  processing  conditions 
the  chemical  pathway  for  the  formation  of  silicon  nitride  is  elemental  silicon,  undergoing 
sublimation  and  a  subsequent  gas  phase  reaction  with  nitrogen: 

SSi,  +  2N2,g^{Si3N4)s  (1) 

where  the  subscripts  g  and  s  indicate  the  gaseous  and  solid  nature  of  the  reactants  and 
the  product.  This  process  is  essentially  similar  to  a  chemical  vapor  deposition  process  with 
product  formation  on  previously  formed  silicon  nitride.  The  free  surface  of  silicon  nitride 
serves  as  a  reaction  catalyst.  Although  the  complete  reaction  may  include  intermediate 
steps  such  as  formation  of  silicon  monoxide  vapor  which  then  acts  as  a  silicon  carrier,  here 
these  details  are  not  of  interest  but  rather  a  single  all-encompassing  composite  reaction  is 
developed.  The  formation  of  a  reaction  product  on  previously  existing  silicon  nitride  gives 
rise  to  a  growing  product  front  which  is  in  competition  for  space  with  the  silicon  particles, 
in  the  process  of  receding  by  evaporation.  The  balance  between  these  two  processes  is  one  of 
the  phenomena  which  most  strongly  affects  the  reaction  rate  and  morphology  of  the  unique 
porous  microstructure. 

Silicon  transport  by  solid  state  diffusion  has  been  suggested  by  some  investigators  [6] 
however  this  is  kinetically  unfavorable.  Solid  state  diffusion  is  too  slow  to  account  for  the 
observed  “fast”  reaction  rates  [7,8].  In  addition,  it  has  been  pointed  out  that  the  silicon  vapor 
pressure  of  10~^  atmospheres  at  typical  processing  temperatures,  is  more  than  adequate  to 
sustain  the  observed  nitridation  rates  [9]. 

The  assumption  of  fast  nitrogen  vapor  transport  has  been  verified  experimentally  [1].  For 
1  atmosphere  N2  and  samples  two  millimeters  thick,  complete  nitridation  has  been  observed 
at  the  center  of  samples  of  average  thickness  in  the  millimeter  range  for  temperatures  from 
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1150  to  1350°C  over  typical  reaction  periods  of  5  to  60  minutes  [10].  Therefore,  in  the  present 
simulation,  a  constant  partial  pressure  of  nitrogen  is  used,  and  the  silicon  vapor  is  the  rate 
limiting  medium,  adjusting  to  new  conditions  instantaneously. 

Popper  and  Ruddlesden  [11]  recognized  that  the  exothermic  nature  of  the  reaction  could 
cause  melting  of  the  unreacted  silicon  unless  the  heating  rate  was  carefully  controlled.  How¬ 
ever,  such  melting  effects  are  considered  inessential  and  have  not  been  incorporated  in  the 
present  model.  Isothermal  conditions  are  assumed  for  simplicity. 

Lastly,  acknowledging  that  the  silicon  nitride  grains  exhibit  anisotropic  growth,  we  have 
allowed  for  preferential  growth  directions.  The  degree  of  anisotropy  is  dependent  upon  which 
of  two  silicon  nitride  phases  dominate;  a  or  While  the  precise  control  of  phase  content 
is  currently  a  topic  of  debate,  it  has  been  observed  that  the  a  phase  is  favored  for  vapor 
phase  reaction  processes  at  lower  temperatures,  and  for  relatively  pure  silicon  powders.  Its 
morphology  has  been  observed  to  range  from  equiaxed  to  whisker- like.  The  7?  phase  has  an 
elongated  morphology,  though  not  whisker-like,  and  is  typically  larger  in  size.  It  is  prominent 
in  materials  processed  with  either  impurities  such  as  iron  or  at  higher  temperatures  where 
liquids  are  prevalent. 

2.2  The  Numerical  Technique 

Cellular  automata  methods  have  been  developed  for  modeling  microscopic  physical  pro¬ 
cesses  to  study  fluid  dynamics  [12,13],  solid  state  diffusion  [14]  and  critical  points  in  phase 
transitions  [15].  Recently,  this  technique  has  also  been  used  to  address  the  issues  in  the  com¬ 
plex  coupled  problems  of  reaction/diffusion  associated  with  materials  processing.  Dab  and 
Boon  [16]  have  presented  results  for  a  more  sophisticated  microscopic  process  than  the  usual 
simple  collisions  in  the  previous  models,  by  introducing  the  existence  of  chemical  reactions  in 
the  medium.  The  model  of  an  autocatalytic  reaction  incorporating  diffusion,  displayed  essen¬ 
tial  features  in  qualitative  agreement  with  systems  involving  non-discrete  reaction- diffusion 
as  described  by  analytical  solutions  of  such  processes.  More  recently,  Srolovitz  and  co-workers 
have  demonstrated  good  agreement  between  experiment  and  three  dimensional  Monte  Carlo 
simulations  of  grain  coarsening  [17]. 

In  the  present  study  the  modeling  of  the  reaction  has  been  performed  in  the  spirit  of  the 
cellular  automata  (CA)  approach,  where  the  simulation  volume  is  divided  into  cells,  each  of 
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which  is  allowed  to  exist  in  one  of  few  discrete  states.  The  state  of  a  given  cell  changes  in 
steps,  as  determined  by  its  current  state  and  that  of  its  nearest  neighbors.  However,  unlike 
the  classical  CA  approach,  where  evolution  of  states  is  deterministic,  here  the  cell  in  any 
given  state  has  more  than  one  possible  transformation  alternative. 

To  choose  between  the  several  possible  outcomes  in  the  reaction,  a  random  number 
generator  was  employed.  Thus,  the  approach  used  here  has  features  common  to  both  CA 
and  Monte  Carlo  algorithms.  However,  deterministic  CA  evolutions  may  be  thought  of  as 
limiting  cases  of  the  present  stochastic  simulation,  when  probabilities  of  all  but  one  outcome 
are  set  to  zero.  This  is  why  we  associate  the  present  simulations  with  the  CA  approach. 

The  reaction  model  that  has  been  used  is  a  2-D  cut  out  of  a  3-D  dense  random  packing  of 
spheres  in  space,  as  will  be  expanded  on  below  in  Section  2.3.  The  choice  of  resolution  of  the 
mesh  was  based  on  maintaining  equiaxed  particle  shapes,  therefore,  the  cell  dimension  was  a 
small  fraction  of  the  initial  silicon  particle  size.  The  total  mesh  dimension  was  limited  only  by 
computational  resources.  At  maximum  resolution,  we  have  obtained  fully  reacted  structures 
spanning  1000  by  1000  cells  or  approximately  130  particle  diameters  or  equivalently  a  volume 
of  2  X  10®  particles.  In  order  to  expedite  generation  of  larger  number  of  structures,  a  mesh 
of  an  edge  length  of  64  cells  was  chosen  for  most  simulations.  This  insures  the  packing  of 
enough  silicon  particles  into  the  mesh  (about  100)  and  provides  an  acceptable  representative 
volume  element  (RVE)  for  a  relevant  microstructure. 

2.3  Motivation  and  Limitations  of  the  Model 

As  already  stated,  the  motivation  for  developing  the  present  numerical  model  was  to  gen¬ 
erate  microstructures  for  subsequent  fracture  studies,  rather  than  to  develop  a  fully  detailed 
study  of  the  reaction  kinetics.  The  limitations  of  the  model,  primarily  its  two  dimensional 
nature  and  its  fixed  resolution,  restrict  our  conclusions  to  indications  of  trends.  Neverthe¬ 
less,  it  has  been  possible  to  obtain  meaningful  trends  of  fracture  toughness  that  could  be 
identified  with  specific  processing  parameters.  In  addition,  it  appears  possible  to  obtain 
some  insight  into  such  complex  phenomena  as  the  spatial  competition  between  the  growing 
silicon  nitride  particles  and  the  silicon  particles  receding  by  evaporation  of  silicon,  and  the 
elfects  these  processes  have  on  the  fracture  toughness  of  the  material. 

In  the  present  model  any  internal  stresses  arising  from  such  effects  as  densification  through 


6 


sintering  or  volume  changes  arising  from  the  reaction,  other  than  from  the  stoichiometry  of 
the  chemical  process,  are  not  specifically  considered.  Incorporation  of  this  effect  would 
have  required  a  substantial  expansion  of  the  model  to  include  simultaneous  considerations 
of  mechanical  equilibrium  and  elastic  strain  energies.  We  compensate  for  this  omission  by 
assuming  certain  parameters  governing  the  strength  of  interfaces  in  our  fracture  simulation 

[4]. 
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3.  THE  MODEL 


3.1  Overall  Outline 

The  following  features  of  the  model  will  be  described  operationally:  the  packing  of  the 
unreacted  green  body;  the  placement  of  initial  silicon  nitride  nuclei  which  become  single 
grains  in  the  final  microstructure;  and  finally  the  chemical  conversion  into  reaction  bonded 
silicon  nitride  (RBSN). 

3.2  Packing  of  Unreacted  Silicon  Particles 

To  avoid  artificial  clustering  phenomena  and  excessive  porosity  arising  from  the  two  di¬ 
mensional  nature  of  our  model,  the  initial  unreacted  silicon  particles  were  actually  modeled 
as  a  three  dimensional  random  close  packed  structure  of  uniform  sized  spheres.  This  ap¬ 
proximation  is  well  suited  to  the  actual  spherical  Si  particles  used  in  the  RBSN  produced  by 
Haggerty  and  co-workers  [1,2]  and  for  its  qualitative  similarity  to  the  usual  green  body  pack¬ 
ings  (i.e.  the  existence  of  expected  naturally  occurring  packing  defects  and  lack  of  long-range 
periodicity  [18])  with  its  similar  packing  density,^  (62%). 

The  particle  coordinate  data  was  generated  by  a  molecular  dynamics  (MD)  technique 
using  periodic  boundary  conditions.  Many  methods  can  be  utilized  to  generate  dense  random 
packing  of  spheres,  but  the  MD  technique  [19]  is  efficient  and  provides  packings  which  have 
the  spatial  distribution  of  the  random  close  packed  structure  as  demonstrated  by  the  radial 
pair  distribution  function  of  particle  centers  in  Fig.  2. 

This  3-D  structure  of  dense  random  packed  spheres  was  subsequently  sectioned  by  a 
geometrical  plane  to  obtain  a  two  dimensional  arrangement  of  circles  of  a  well  known  size 
distribution.  This  information  was  then  discretized  using  on  a  hexagonal  planar  mesh. 

Figure  3  shows  a  typical  planar  section  through  the  3-D  initial  dense  random  packing 
of  spheres  which  was  taken  for  the  2-D  model.  Figure  4  shows  schematically  how  the  2-D 
plane  arrangement  of  material  in  circular  disc  form  and  surrounding  voids  are  discretized 

^  Dense  random  structures  have  the  most  efficient  packing  of  uniform  sized  particles  without  long  range 
periodicity. 
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into  reacting  cells  for  the  reaction  model. 


3.3  Nucleation  of  Silicon  Nitride 

Once  the  geometry  of  the  sectioned  particle  field  is  mapped  on  the  2-D  mesh  hexagonal, 
the  cells  of  the  mesh  that  represent  surface  sites  of  Si  particles  are  located  (i.e.  Si  cells 
with  at  least  one  neighboring  empty  site).  A  given  number  of  these  surface  sites  were  then 
randomly  chosen  to  be  the  initial  Si3N4  nuclei  and  were  replaced  by  S'i3A^4.  The  number  of 
such  nuclei  was  small,  typically  1-10  per  particle,  (i.e.  100-1000  nuclei  compared  to  typically 
10,000  total  cells). 

In  some  instances,  a  given  proportion  of  these  nuclei  were  assigned  to  have  preferential 
growth  directions  to  simulate  anisotropy  of  grain  growth.  Growth  directions  were  chosen  ran¬ 
domly  from  among  three  discrete  directions  corresponding  to  the  three-principal  symmetry 
axes  directions  of  the  mesh. 

3.4  Conversion  Reaction  to  Si^N^ 

3.4-1  Basic  components 

Given  the  green  structure  of  Si  particles,  the  probabilities  were  computed  for  each  of 
the  following  possible  events  at  each  individual  surface  cell  site:  (a)  the  evaporation  and 
condensation  of  silicon;  (b)  the  reaction  and  growth  of  silicon  nitride;  (c)  grain  coarsening 
and  surface  diffusion. 

3.4.2  The  condensation  and  evaporation  of  silicon  vapor  to  and  from  free  silicon  surfaces 

For  every  site  that  represents  Si,  each  nearest  neighbor  was  checked  for  an  empty  cell.  If 
one  existed,  the  Si  site  was  labeled  as  a  surface  Si  site,  and  was  a  possible  candidate  for  evap¬ 
oration.  The  incremental  probability  for  evaporation,  ^F’evap,  was  taken  to  be  proportional 
to  the  increment  in  time,  St,  chosen  for  all  processes,  i.e. 

(5Pevap  =  E-St,  (2) 
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where  is  a  kinetical  rate  constant  discussed  further  below.  If  evaporation  was  deemed 
to  occur,  the  cell  was  transformed  into  an  empty  cell  and  the  silicon  vapor  pressure  was 
incremented. 

Condensation  was  taken  to  be  the  reverse  of  evaporation  for  which  an  empty  cell  located 
next  to  a  Si  cell  was  a  candidate  for  condensation,  where  the  probability  for  the  condensation 
transformation  was  taken  to  be  proportional  to  the  silicon  vapor  pressure.^ 

3.4-3  The  conversion  reaction  and  growth  of  Si3N4  from  previously  grown  free  Si3N4  surfaces 

In  a  similar  manner  as  for  the  evaporation  and  condensation  processes  presented  above, 
empty  cells  neighboring  on  Si3N4,  were  allowed  the  possibility  of  undergoing  a  transformation 
to  5^3 AI4,  thereby  extending  the  growth  front.  The  incremental  probability,  ^Preactiom  of 
growth  of  Si3N4,  was  taken  to  be  proportional  to  the  partial  pressure  Psiy  of  silicon  vapor 
(the  rate  limiting  species),  a  rate  constant  R,  and  the  time  step,  8t: 

reaction  ~  P  ’  Psiv  '  8t.  (3) 

Once  growth  occurred,  the  vapor  phase  was  reduced  in  a  stoichmetrically  appropriate 
way.^  It  is  important  to  note  that  this  generic  growth  law  does  not  specify  the  mechanism  of 
growth.  Only  that  previously  grown  Si3N4  is  allowed  to  continue  to  grow  further,  and  that 
the  growth  rate  is  proportional  to  the  concentration  of  the  critical  vapor  species  of  silicon,  as 
well  as  the  reaction  constant  R.  It  is  assumed  that  N2  is  present  in  abundance  as  an  infinite 
reservoir,  requiring  no  separate  N2  balance. 

3.4.4  Grain  boundaries  and  surface  tension 

Surface  tension  and  “grain  boundary  tension”  were  incorporated  into  the  model  in  order 
to  give  more  realistic  appearing  structures.  Although  sintering  processes  are  thought  not 
to  be  very  influential  in  the  low  temperature  processing  regime  (T  <  1350°C),  micrographs 
(both  TEM  and  SEM)  give  the  impression  of  relatively  smooth  grain  boundaries  and  pore 

^The  amounts  of  silicon  evaporation  and  condensation  are  implicitly  a  function  of  the  amount  of  silicon 
free  surface  because  each  individual  surface  particle  is  given  a  chance  to  evaporate  or  adsorb  silicon  vapor. 

^The  volumetric  expansion  of  Si  reactant  in  the  reaction  bonding  process  with  N  gas  to  final  Si3N4  is 
1.22.  Therefore  in  order  to  conserve  Si  mass,  to  produce  1  volumetric  unit  of  523/^4  product,  0.82  fractional 
units  of  one  solid  cell  of  Si  is  effectively  consumed,  from  the  vapor  phase. 


10 


surfaces.  Therefore,  surface  tension  and  boundary  tension  parameters  were  introduced  for  the 
silicon  nitride  material.  For  grain  boundary  motion,  the  silicon  nitride  cells  with  neighboring 
cells  of  silicon  nitride  of  another  grain  orientation  were  given  a  finite  probability  to  transform 
to  the  orientation  of  the  other  grain.  In  the  case  of  the  free  surface,  a  surface  cell  was  allowed 
the  possibility  of  translating  to  a  neighboring  empty  surface  site  that  is  adjacent  to  the  same 
grain  of  silicon  nitride.  The  specific  details  are  presented  in  Appendix  I. 

3.4-5  Implementation  of  the  reaction  model 

In  practice,  the  simulation  began  by  sectioning  of  the  3-D  particle  packing  and  mapping 
the  resulting  cross-section  onto  the  2-D  mesh  followed  by  the  placing  of  the  523^4  nuclei 
as  outlined  above.  Once  this  was  accomplished,  the  global  variables,  time  t,  and  silicon 
vapor  pressure,  Psiv  were  initialized  to  t  =  0  and  =  0.  After  this  initialization,  time  was 
incremented  by  St.  The  choice  of  St  was  made  by  computing  incremental  probabilities  for  all 
processes  and  by  making  St  to  be  a  small  fraction  of  the  inverse  rate  constant  of  the  fastest 
process.  In  this  manner,  all  processes  were  scaled  consistently,  so  that  the  transformation 
of  a  particular  cell  was  maintained  to  be  exclusively  a  function  of  the  current  state  of  its 
surrounding  neighbors.  With  this  choice  of  St,  for  each  cell,  all  incremental  transformation 
probabilities  were  computed  and  certain  events  were  chosen  to  take  place,  all  based  upon 
the  use  of  a  random  number  generator.  In  the  cases  of  evaporation,  condensation  and 
reaction  of  Si3N4,  the  mass  balance  of  Si  was  accounted  for  by  appropriately  incrementing 
and  decrementing  the  parameter  Ps,„.  As  already  stated,  the  transport  of  N2  was  assumed  to 
be  easy  in  comparison  to  Si,  so  that  no  mass  balance  was  kept  for  it.  Moreover,  isothermal 
conditions  were  assumed  requiring  no  adjustments  of  the  rate  constants  in  the  course  of  the 
simulations. 
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4.  RESULTS 


4.1  Overview  of  Results 

Many  features  of  the  reaction  kinetics  process  are  mirrored  in  these  simulations.  For 
example,  the  basic  trends  of  the  extent  of  the  reaction  with  time  and  the  pore  to  particle 
size  ratio  are  in  quantitative  agreement  with  experiments  as  we  will  demonstrate  below. 
Moreover,  the  effect  of  the  initial  Si  particle  positions  on  the  final  pore  positions  is  also  in 
qualitative  agreement  with  experimental  observations. 

Making  due  allowances  for  some  noted  assumptions,  we  utilize  the  fact  that  the  simu¬ 
lations  presented  here  have  some  connection  with  the  real  process,  and  report  observations 
of  the  effects  of  the  reaction  on  microstructure.  Throughout  the  discussion  below,  attention 
will  be  paid  to  microstructural  features  as  they  relate  to  the  various  reaction  parameters,  in 
particular  grain  and  pore  size  and  shape,  as  well  as  the  effect  of  packing  defects  on  the  final 
microstructure. 

4.2  Comparison  with  Experimental  Observations 

In  Fig.  5,  the  comparison  of  the  normalized  number  density  distribution  of  pore  sizes  (in 
units  of  particle  size)  for  the  simulated  and  actual  microstructures  reveals  comparable  average 
dimensions  and  qualitative  agreement  between  the  size  distributions.  Other  comparisons  of 
pore  aspect  ratio  and  pore  circumference  distributions  between  the  model  and  experiments 
is  given  in  Fig.  6.  The  ratio  of  the  average  grain  size,  D,  to  the  initial  Si  particle  size,  d, 
for  initial  density  of  nuclei,  N  =  14.4,  per  Si  particle,  in  a  2-D  cross  section  is  0.45  in  the 
model,  compared  to  the  expected  ratio  of  0.439  as  presented  in  Appendix  II.  The  ^24^4  grain 
and  initial  Si  particle  sizes  were  determined  experimentally  with  a  combination  of  electron 
microscopy  and  x-ray  techniques  [20]. 

Additional  agreement  between  the  model  and  actual  reaction  kinetics  behavior  can  be 
found  in  the  time  dependence  of  the  extent  of  the  reaction.  The  simulated  extent  of  the 
reaction  follows  the  Johnson-Mehl-Avrami  form  of: 

X  =  1  —  exp{aNP)  (4) 
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where  X  is  the  total  fraction  transformed,  N  is  the  constant  density  of  nuclei,  a  includes 
the  growth  rate  constant  and  geometric  factors,  t  is  time,  and  n  is  the  time  exponent. 
Figure  7a  shows  the  simulated  dependence  of  X  on  t,  which  compares  well  with  the  form  of 
the  Johnson-Mehl-Avrami  model,  giving  an  exponent  of  ra  =  1.5.  Figure  7b  shows  thermo- 
gravimetric  data  by  Sheldon  and  Haggerty  [1],  given  here  for  comparison  with  the  simulation 
data.  A  key  difference  between  the  simulation  and  experiments  is  in  the  nucleation  regime 
with  its  varying  waiting  times,  followed  by  the  sigmoidal  trend  of  the  JMA  model  in  the 
actual  experiments,  which  were  discussed  by  Sheldon  and  Haggerty  [10].  In  a  possible  further 
refinement  of  the  model  such  incubation  times  representing  nucleation  barriers  could  have 
been  readily  incorporated  by  effectively  re-scaling  the  origin  of  time  in  a  stochastic  manner. 
Here  this  was  not  of  interest,  since  the  principal  goal  was  not  modeling  the  reaction  in  all 
its  detail  but  only  to  generate  non-arbitrary  microstructures.  The  fact  that  this  simple 
model  follows  this  form  and  produces  realistic  microstructures  reinforces  the  notion  that  the 
model  incorporates  correct  key  features  of  the  reaction  bonding  process  -  albeit,  without  a 
nucleation  barrier. 

In  the  following  sections  the  trends  of  pertinent  microstructural  features  obtained  by  the 
simulation  are  presented  with  dimensionless  ratios  of  the  reaction  parameters,  followed  by 
a  discussion  of  the  mechanisms  which  appear  to  control  the  formation  of  large  pores  in  the 
final  microstructure. 

4.3  Reaction  of  Si^N^  vs  Evaporation  of  Si  Rates 

The  competition  between  the  reaction  rate  of  Si  and  N2  and  the  rate  of  evaporation  of  Si 
on  active  surfaces  is  the  most  relevant  consideration  in  the  topological  competition  of  growing 
Si3N4  and  the  receding,  evaporating  Si  particles.  A  large  ratio  of  reaction  rate  constant,  R, 
to  evaporation  constant,  E,  subsequently  referred  to  as  R/E,  indicates  immediate  reaction 
to  form  silicon  nitride  and  a  high  degree  of  accentuation  of  the  kinetics  compared  with  purely 
geometrical  impingement  of  the  silicon  nitride  on  the  more  slowly  receding  silicon. 

For  this  reason  of  increased  impingement  of  the  reaction  product  on  the  Si  particles 
at  high  reaction  rates,  there  is  some  correlation  between  the  positions  of  the  initial  silicon 
particles  and  final  pores  of  sizes  comparable  to  the  initial  particle  size.  This  is  illustrated  in 
Fig.  8a  and  8b  for  two  extreme  R/ E  ratios  of  0.3  and  300.  Figure  8c  gives  for  comparison 
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the  initial  position  of  the  Si  particles.  The  above  mentioned  effect  of  replacement  of  particles 
with  pores  is  evident  in  the  case  with  the  high  R/E  ratio  of  300. 

Large  pores  in  the  initial  packing  are  observed  to  be  the  locations  of  large  final  pores  for 
both  high  and  low  reaction  rates.  The  packing  of  the  initial  particles  in  Figs.  8a  and  8b, 
unlike  the  other  simulations,  were  performed  in  two  dimensions  instead  of  in  three  dimensions 
as  previously  described.  This  was  done  to  illustrate  the  agglomeration  effects  of  packing  of 
circular  shapes  of  equal  diameter  in  2-D,  and  the  associated  larger  packing  defects  which 
arise.  However,  the  high  reaction  rate  microstructure  to  some  extent  has  “healed”  the  largest 
pores  associated  with  the  large  packing  defects.  This  is  because  the  Si3N4  growth  has  been 
stifled  by  its  impingement  with  neighboring  Si  particles  everywhere  except  where  the  growth 
front  had  ample  space  to  grow,  i.e.,  at  the  large  packing  defects.  Thus,  conditions  with  large 
R/E  provide  for  high  growth  rates  for  those  nuclei  which  are  growing  unimpeded  in  the  space 
of  the  large  voids.  This  effect  is  shown  clearly  in  the  shrinking  of  the  large  pore  in  the  upper 
right  quadrant  of  Fig.  8b  and  in  Figs.  9a  and  9b  which  show  representative  microstructures 
of  a  larger  scale  for  the  two  cases  shown  in  Figs.  8a  and  8b,  respectively. 

Finally,  the  amount  of  residual  trapped  silicon  increases  for  high  R/E  ratio  conditions 
due  to  increased  impingement  of  Si3N4  with  the  attendant  encapsulation  of  the  evaporating 
Si,  and  the  reduction  of  the  Si  surface  available  for  evaporation,  as  in  the  case  of  Fig.  8b 
for  the  R/E  ratio  of  300.  Note  that  the  unreacted  Si  is  located  in  the  most  densely  packed 
region  of  the  initial  particle  packing,  where  impingement  is  expected  to  be  the  greatest. 
Conversely,  a  low  value  of  the  R/E  ratio  implies  that  ample  silicon  vapor  is  created  and 
impingement  is  reduced.  Figures  10a  and  10b  show  the  separate  process  kinetics  plots  for 
R/E  of  0.3  and  300.  Clearly,  a  decrease  in  the  reaction  rate  R  by  1000  at  constant  E,  which 
spans  the  transition  from  impinging  to  non- impinging  reactant /product  fronts,  results  in 
radical  differences.  For  the  case  of  the  R/E  ratio  of  300  where  Si  evaporation  is  followed 
by  immediate  consumption  of  the  vapor  by  the  silicon  vapor  starved  reaction,  the  time  for 
completion  of  the  reaction  is,  nevertheless,  shortened  to  only  1/4  of  the  time  for  the  case  of 
R/E  ratio  of  0.3. 
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4.4  Evaporation  vs  Condensation  Rates 


A  high  evaporation  to  condensation  rate  ratio  E  jC  produces  higher  vapor  pressures  and 
a  proportionately  higher  reaction  rate.  This  effect  is  quite  strong  as  can  be  seen  in  Figs. 
11a  and  11b  for  the  EjC  ratios  of  0.1  and  10  respectively.  For  a  100  fold  increase  in  the 
EjC  ratio,  with  all  other  parameters  held  constant,  the  time  for  completion  of  the  reaction 
is  shorter  by  a  factor  of  300.  The  additional  amplification  by  a  factor  of  3  of  this  effect 
appears  to  be  due  to  reduced  impingement  of  growing  Si:iN4  interfaces  with  the  Si  particles. 
Impingement  can  be  reduced  if  the  Si  reactant  can  be  effectively  “stored”  in  the  mobile  vapor 
phase.  Figure  12  shows  the  unreacted  Si  in  the  microstructure  resulting  from  the  increased 
impingement  for  a  low  evaporation  rate. 

4.5  Density  of  Si^N^  Nuclei  on  Particles 

The  number,  N,  of  initial  reactive  nuclei  on  particles,  influences  the  reaction  rate,  the 
amount  of  unreacted  Si  and  both  the  scale  and  the  homogeneity  of  the  final  microstructure. 
At  the  initial  stages  of  the  reaction,  a  large  number  of  nuclei  per  particle  will  increase 
the  reaction  rate  due  to  the  larger  number  of  reaction  sites.  However,  once  a  comparable 
amount  of  Si^N^  is  generated,  the  growth  rates  are  equivalent.  This  effect  is  shown  in 
Fig.  13  in  terms  of  reaction  history,  microstructures,  and  distribution  of  grain  and  pore  size 
distributions,  for  two  different  cases  of  the  initial  number  of  Si^N^  nuclei  per  particle  for  the 
same  RjE  and  EjC  ratios  of  3  and  1,  respectively.  In  Fig.  13a  the  reaction  history  and  the 
resulting  distributions  of  grain  and  pore  sizes  together  with  the  associated  microstructure 
are  shown  for  a  case  of  high  density  (10)  of  nuclei  per  particles  while  Fig.  13b  show  the  same 
corresponding  information  for  a  low  density  of  nuclei.  Extremely  high  density  of  nuclei  for 
which  the  Si  evaporation  is  blocked,  stops  the  reaction  prematurely. 

If  the  density  of  nuclei,  N,  is  low  enough,  fewer  than  10  per  particle  in  our  observa¬ 
tions,  and  the  nuclei  are  distributed  randomly,  there  is  a  statistically  relevant  probability  for 
forming  green  structure  regions  with  zero  nuclei,  resulting  in  large  voids  in  the  final  reacted 
structure.  Such  a  case  is  shown  in  Fig.  14  and  Fig.  13b,  with  large,  interconnected  pores 
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having  tortuous  borders.  Conversely,  a  large  density  of  nuclei,  N,  reduces  this  probability 
and  distributes  the  reaction  product  more  uniformly  as  is  clear  from  Fig.  13a.  Note  that 
the  number  of  grains  or  nuclei  per  particle  in  real  microstructures  is  observed  to  be  approx¬ 
imately  7  in  a  2-D  cross  section,  which  is  close  to  the  critical  density  for  formation  of  large 
scale  voids. 

Finally,  with  fewer  nuclei,  the  final  average  grain  size  is  larger,  because  each  nucleus 
at  the  initial  stage  represents  a  grain  in  the  final  stage,  with  the  total  amount  of  reacted 
material  divided  among  those  fewer  grains.  Quantitatively,  the  pore  size  scales  with  grain 
size,  which  becomes  evident  by  the  average  grain  to  pore  size  ratio  which  remains  constant 
(roughly  at  a  level  of  2)  over  a  range  of  density  of  nuclei,  N,  from  1.0  to  10  per  particle. 

4.6  Packing  Defects 

Shapes  of  larger  pores  are  observed  to  be  related  to  three  types  of  defects  in  the  packing 
of  unreacted  particles.  Some  of  these  defect  types  are  shown  in  Fig.  15,  and  include:  (1) 
regions  of  high  density  of  packing,  where  there  is  highly  restricted  growth,  as  shown  in  region 
A  of  Fig.  15;  (2)  large  scale  pores,  on  the  order  of  multiple  particle  diameters  as  shown  in 
region  B  of  the  figure,  where  an  initially  low  density  of  particles  results  in  a  large  final  pore; 
or  (3)  low  density  of  nuclei  which  was  discussed  above.  All  three  types  of  defects  have  been 
observed  in  the  simulations  to  produce  large  pores  in  the  final  microstructure.  The  first  and 
second  cases  are  shown  to  coexist  in  Fig.  15,  which  shows  the  resulting  microstructure  with 
and  without  the  particle  overlay. 

There  is  reason  to  believe  that  the  large  scale  packing  defects  in  the  green  structure  such 
as  cases  1  and  2  above  exist  in  real  structures.  Large  samples  of  uniform  size  powders  are 
likely  to  be  aggregates  of  dense  and  highly  ordered  crystalline  domains,  with  a  lower  density 
of  randomly  packed  particles  in  between  [21,22].  A  key  difference  between  the  ideal  random 
close  packing  procedure  and  the  real  pressed  green  compacts  is  that  the  latter  experience 
interparticle  friction  and  associated  concentrations  of  pressure  resulting  in  possible  localized 
densifi  cation. 
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5.  DISCUSSION 


The  simulation  presented  here  primarily  provides  an  objective  means  for  generating  real¬ 
istic  appearing  microstructures  of  reaction  bonded  porous  silicon  nitride.  The  fact  that  the 
grain  and  pore  sizes  and  shapes  appear  to  closely  resemble  those  of  real  microstructures  and 
that  the  simulation  mirrors  other  features  of  the  reaction  process,  places  more  credence  on 
this  technique.  Moreover,  the  physical  basis  of  the  model  implies  that  these  microstructures 
can  be  altered  in  an  objective  manner  by  varying  parameters  which  are  identifiable  with 
physical  quantities  of  the  reaction.  In  addition,  it  can  be  noted  that,  it  is  possible  to  obtain 
some  insight  into  such  complex  phenomena  as  the  spatial  competition  between  the  grow¬ 
ing  silicon  nitride  and  the  receding,  evaporating  silicon.  The  observed  reaction  mechanisms 
depend  not  only  on  elementary  reaction  rates,  but  on  these  structural  effects. 

The  most  prominent  parameters  in  the  model  which  govern  the  microstructural  homo¬ 
geneity,  pore  and  grain  shape  and  size  are  the  reaction  to  evaporation  ratio  Rj E,  the  evapo¬ 
ration  to  condensation  ratio  EjC^  and  the  initial  density  of  silicon  nitride  nuclei  on  particle 
surfaces.  All  of  these  parameters  directly  affect  the  degree  of  product/reactant  impingement, 
the  time  for  completion  of  the  reaction  and  the  homogeneity  of  the  microstructure.  Specif¬ 
ically,  for  the  cases  of  low  density  of  nuclei,  high  Rj E  ratio  and  low  EjC  ratio,  significant 
product  impingement  is  found  to  occur.  In  this  regime,  the  following  effects  are  observed: 
reduced  overall  reaction  rates,  a  high  correlation  between  initial  Si  particle  positions  and 
final  pore  positions,  and  a  tendency  for  faster  product  growth  in  regions  of  lower  packing 
density,  i.e.,  “healing”  of  packing  defects. 

Several  severe  assumptions  limit  the  present  method  of  simulation  of  the  complex  reaction 
processes  to  the  creation  of  schematic  processing  maps.  Foremost,  is  a  lack  of  specific 
functional  forms  for  reaction  rates.  Therefore,  we  have  tailored  the  analysis  in  terms  of 
non-dimensional  ratios  of  specific  reaction  constants  and  limited  the  results  to  trends.  Also, 
the  limit  of  resolution  of  the  mesh  and  the  two  dimensional  nature  of  the  model  place  some 
degree  of  artificial  constraints  on  the  simulated  processes.  The  resolution  of  the  mesh  is  less 
problematic  if  the  continuum  processes  can  be  thought  of  occurring  smoothly  during  the 
discrete  time  increments  and  if  the  typical  microstructural  feature  is  sufficiently  larger  than 
the  mesh  resolution.  The  two  dimensional  nature  of  the  simulation  has  been  compensated 
for  partly  by  using  a  2-D  section  of  a  true  3-D  packing  of  particles  for  the  starting  condition 
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and  by  appropriately  adjusting  the  density  of  S'i3A^4  nuclei  to  conserve  the  grain  size  to 
Si  particle  size.  Since  the  microstructure  is  certain  to  evolve  in  all  three  dimensions  in  a 
statistically  similar  way,  we  have  compensated  for  the  lack  of  three  dimensionality  in  an 
effective  manner  in  the  simulation  of  crack  propagation  through  the  microstructure  by  the 
introduction  of  an  effective  3-D  element  based  on  a  specific  grouping  of  near  neighbors  of 
microstructure  elements.  This  we  describe  in  the  companion  communication  [4]. 

Two  situations  are  too  complex  to  simulate  to  justify  the  relatively  small  effects  which 
would  be  predicted  to  occur.  The  first  is  the  isolation  of  large  regions  to  either  Si  vapor  or  Ni 
gas  reactants.  The  approximation  of  small  Si  vapor  gradients  is  valid  when  the  evaporation 
of  the  Si  particles  occurs  homogeneously  throughout  the  reaction.  However,  if  structurally  a 
reaction  product  (or  unreacted  Si)  forms  a  barrier  around  a  volume  larger  than  one  particle, 
which  can  be  accounted  for  by  the  simulation,  the  reaction  within  the  region  will  naturally 
be  choked  off.  Because  in  the  typical  reaction  case,  nearly  complete  nitridation  is  reached 
[10],  the  likelihood  of  the  occurrence  of  many  reactant  starved  regions  is  minimal.  A  second 
situation  of  considerable  complexity  which  we  have  not  considered  is  the  development  of 
reaction  induced  internal  stresses  of  misfitting  material.  This  requires  a  simultaneous  solu¬ 
tion  of  mechanical  equilibrium  and  introduces  interactions  such  as  the  effect  of  pressure  on 
reaction  rates.  In  spite  of  the  additional  insight  this  would  have  provided,  the  associated 
complexities  outweigh  the  advantages.  In  the  companion  fracture  simulation  we  have  com¬ 
pensated,  in  part,  for  the  absent  effects  of  development  of  internal  stresses  by  introducing 
material  strength  and  stiffness  variability. 

With  these  observations  and  the  ability  to  generate  unique  microstructures  which  would 
be  difficult  (if  not  impossible)  to  produce  or  precisely  vary  in  the  laboratory,  one  can  then 
proceed  to  sample  the  fracture  properties  of  these  microstructures.  In  this  way,  not  only  can 
trends  of  fracture  toughness  be  associated  with  processing  variables,  but  also  a  mechanistic 
understanding  of  the  reasons  for  the  observed  trends  can  be  brought  to  the  material  design 
problem  to  introduce  pre-selected  microstructural  design  criteria  for  subsequent  material 
processing  iterations.  The  results  of  the  companion  fracture  paths  study,  however,  show 
that  the  degrees  of  effective  freedom  are  not  larger  in  manipulating  material  toughness.  [4] 
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6.  CONCLUSIONS 


The  kinetics  simulation  presented  here  captures  many  of  the  salient  features  of  the  reac¬ 
tion  bonding  process  in  terms  of  considerations  for  determining  fracture  toughness.  Primar¬ 
ily,  the  generated  microstructures  resemble  actual  reaction  bonded  silicon  nitride  microstruc¬ 
tures,  in  terms  of  the  grain  and  pore  to  particle  size  ratios.  The  observed  microstructural 
evolution  mechanisms  depend  not  only  on  elementary  reaction  rates,  but  also  on  such  com¬ 
plex  spatial,  structural  phenomena  such  as  the  competition  between  the  growing  SisN^  and 
the  receding,  evaporating  Si  which  are  summarized  below. 

•  Key  dimensionless  reaction  variables  control  the  competition  between  the  growth  of 
the  product  and  the  recession  of  the  reactant,  specifically  the  reaction  to  evaporation 
ratio,  the  evaporation  to  condensation  ratio,  and  the  number  of  initial  S'z3A^4  nuclei 
per  Si  particle. 

•  The  final  pore  positions  are  related  to  the  particle  positions  in  the  unreacted  green 
structure  when  the  product  growth  is  inhibited  due  to  impingement  of  growing  5'z3A^4 
with  the  Si  reactant. 

•  Three  types  of  packing  defects  may  give  rise  to  large  pores  in  the  final  structure: 

1.  Packing  flaws  on  the  scale  of  multiple  particle  dimensions  will  tend  to  remain  as 
large  pores  in  the  final  microstructure.  This  effect  can  be  reduced,  and  defects 
can  be  “healed”  to  some  extent,  if  the  reaction  takes  place  in  a  regime  where  the 
reaction  rate  limit  is  geometrical  and  not  of  a  chemical  kinetics  nature. 

2.  Densely  packed  regions,  which  are  known  to  exist  in  real  structures,  retard  growth 
at  the  early  stages  due  to  strong  impingement. 

3.  Regions  of  low  density  of  nuclei  may  form  statistically  if  the  average  density  of 
nuclei  is  lower  than  some  threshold.  Interstitial  regions  which  have  no  nuclei 
associated  with  them  will  be  shielded  from  the  growing  reaction  product,  and 
could  remain  as  unreacted  inclusions. 
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APPENDIX  I.  SURFACE  TENSION  AND  GRAIN  BOUNDARY 

TENSION  OF  5^3iV4. 

The  excess  interface  energy  AE  is  a  thermodynamic  driving  force  providing  for  the 
smooth  grain  boundaries  appearing  on  the  micrographs.  This  energy  is  essentially  the  source 
of  the  surface  tension  in  the  case  of  a  grain-pore  interface,  and  of  the  grain  boundary  tension 
in  the  case  of  a  grain- grain  interface.  In  both  cases,  substantially  large  AE  will  tend  to 
reduce  the  interface  area  (perimeter  in  the  2-D  case).  It  is  rather  straightforward  to  incor¬ 
porate  such  interface  tension  effects  in  the  present  reaction  model  by  introducing  an  energy 
penalty  for  the  transformations  that  result  in  an  increase  of  the  interface  perimeter. 

More  specifically,  the  diffusion  of  a  grain-grain  interface  is  described  here  as  a  net  result 
of  stochastic  transformations  of  the  silicon  nitride  cells  near  the  boundary  of  a  grain  with 
one  orientation  into  the  silicon  nitride  cells  of  a  neighboring  grain  with  another  orientation. 
If  we  introduce  an  energy  penalty,  per  unit  increment  of  the  grain  boundary  perimeter 
then  the  rates  of  such  mass-exchange  transformations  are  computed  as 

^Pqb  motion  —  C'expl-A^'GfiAnGslf^t,  (AI.l) 

where  Atiqb  is  an  increase  in  the  number  of  the  grain  boundary  segments  associated  with  a 
given  transformation,  C  and  A  are  constants,  and  St  is  a  time  step.  The  minus  sign  in  the 
exponent  will  reduce  probabilities  of  transformations  resulting  in  an  increase  of  the  grain 
boundary  perimeter  {Atigb  >  0);  contrariwise  {Augb  <  0)  will  enhance  the  probabilities  of 
transformation  reducing  the  grain  boundary  perimeter.  Clearly,  such  transformations  will 
tend  to  make  the  grain  boundaries  smoother  in  a  natural  way. 

In  the  case  of  the  free  surface,  a  single  surface  element  is  allowed  a  possibility  of  translating 
to  a  neighboring  empty  site  that  is  adjacent  to  the  same  grain.  Likewise,  the  kinetic  rates 
for  such  translations  are  computed  as 

^■^surface  diffusion  C  exp{— A'I'5An5}(5t,  (AI.2) 

where  An^  is  an  increase  in  the  number  of  the  grain  surface  segments  associated  with  a 
given  translation  and  is  the  energy  per  surface  segment  for  such  an  increase  to  occur. 
By  choosing  appropriate  values  for  C  and  A,  the  rates  given  by  (AI.l)  and  (AI.2)  can  be 
adjusted  relative  to  the  rates  of  the  nitridation  reaction,  evaporation,  and  condensation.  On 
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the  other  hand,  the  two  rates  can  be  made  higher  or  lower  relative  to  each  other  by  choosing 
parameters  ^.nd  appropriately.  Accordingly,  by  varying  these  parameters,  grain 
boundaries  and  surfaces  can  be  made  smoother  or  rougher,  as  needed. 
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Appendix  II.  Density  of  Si^N^  Nuclei 


It  is  desirable  to  relate  the  2-D  equivalent  density  of  nuclei,  {N^ j ,  to  the  3-D 
experimentally  measured  density,  (N^ where  p  stands  for  the  Si  particle,  and  g  stands 
for  the  reacted  Si3N4  grain. 

By  standard  stereology,  when  a  cut  is  taken  through  particles  packed  into  3-D, 


AP  VP 

j^tot  ~  ytot  ’ 

where  tot  indicates  the  total  section  area  and  volume. 

Assuming  that  the  same  is  true  for  the  fully  reacted  structure  we  obtain. 


(AII.l) 


Thus, 


AP  VP 

j^tot  ~  ytot  ■ 


A^_VP_ 
j\p  ~  yp' 


(AII.2) 


(AII.3) 


In  the  2-D  section  through  the  unreacted  Si  green  structure  and  fully  reacted  structure. 


A'’  =  NPaF 


(AII.4) 


AP  =  NPaP 


(AII.5) 


where  Np  and  Np  are  the  number  of  particles  and  grains  respectively,  Ap  and  Ap  are  the 
total  areas  of  the  particles  and  grains  respectively,  and  qP  and  aP  are  the  average  areas  in 
the  sectioned  layer  of  the  particles  and  grains. 

Thus, 

fNP\^^  aP  AP  aPVP  /  a  tt 

(  ATti  )  ~  An  ~  -n  T/n  (AII.6). 


\NpJ  0,9  Ap  opVp  ^  ' 

Then,  using  the  approximation  of  spherical  grains  and  Si  particles,  we  obtain  for  the  ratio 
of  the  maximum  cross-sectional  areas 


yp  \  2/3 


(AII.7) 
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Assuming  further  that  the  ratio  of  the  average  cross-sectional  areas  is  the  same,  it  follows 
that 


aP  _ 

0,3  a3 


Vu^j  “  VsjNs' 


where  and  are  the  volumes  of  Si^N^  grains  and  Si  particles  respectively. 


Combining  equations  (AII.6),  and  (AII.8),  we  obtain  the  result, 


(AII.8) 


vw  ^  [vi^/  J  ^  Iw; 

The  term  (^)  is  the  ratio  of  the  total  volume  of  the  reacted  grains  to  the  total  volume 
of  the  initial  Si  particles,  or  simply  the  reaction  volume  expansion  parameter  1.22. 

For  an  average  grain  size  of  100  nm  and  initial  particle  size  of  250  nm,  we  obtain  a  3-D 
nuclei  density  of  19.1  per  Si  particle  and  using  equation  (AII.9)  we  obtain  an  equivalent  2-D 
density  of  nuclei  of  7.5. 
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Figure  1 
Figure  2 


Figure  3 


Figure  4 


Figure  5 


Figure  6 


Figure  7 


Figure  8 


Representative  microstructure  of  reaction  bonded  silicon  nitride. 

Radial  distribution  function  of  coordinates  for  centers  of  spheres  obtained  from 
a  molecular  dynamics  simulation  of  the  structure  of  dense  random  packing  of 
spheres,  using  the  approach  of  Finney  [19]. 

Distribution  of  dense  randomly  packed  spheres  in  3-D  obtained  by  a  molecular 
dynamics  program,  subsequently  sampled  by  a  plane,  giving  a  2-D  starting 
structure  for  the  RBSN  reaction  simulation. 

The  overlapping  of  the  2-D  arrangement  of  random  circles  by  a  hexagonal  mesh 
defining  the  cells  for  the  cellular  automata  reaction  simulation.  The  material 
inside  the  circles  represents  particles,  black  cells  are  5'z3A^4  nuclei. 

Pore  size  distributions  for  real  and  simulated  microstructures.  The  simulated 
structure  had  reaction  parameters  as  follows:  R/E  =  0.3,  E/C  =  1,N  =  25 
(initial  density  of  nuclei). 

Distributions  of  pore  shape  parameters  for  real,  (a),  and  simulated,  (b),  mi¬ 
crostructures  of  RBSN.  The  simulated  structure  had  reaction  parameters  as 
follows:  R/E  =  0.3,  E/C  =  1,  N  =  7.6. 

(a)  Simulated  history  of  completion  of  the  RBSN  reaction  as  a  function  of  time 
(x  is  the  total  transformed  fraction)  showing  Johnson-Mehl-Avrami  kinetics; 

(b)  experimental  thermo-gravimetric  results  of  Sheldon  and  Haggerty  [1]  for 
the  RBSN  reaction.  The  four  curves  are  for  four  separate  observations. 

Simulated  microstructures  for  low  (R/E  =  0.3):  (a),  high  (R/E  =  300);  (b),  ra¬ 
tios  of  reaction/evaporation  rates,  and  (c)  showing  positions  of  initial  random 
packed  particles.  Packed  particles  were  equal  diameter  discs  packed  in  2-D  to 
illustrate  agglomeration  effects  and  the  tendency  of  the  high  R/E  simulation 
to  “heal”  packing  defects.  (E/C  =  1,  Nuclei  density,  N  =  10.) 
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Figure  9 


Figure  10 


Simulated  RBSN  microstructure:  (a),  for  R/E  =  0.3,  E/C  =  1,  N  =  25;  and 
(b),  for  R/E  =  300,  E/C  =  1,  N  =  25; 

Relative  concentrations  of  various  reaction  species:  (a),  for  low  reac¬ 
tion/evaporation  ratio  (R/E  =  0.3);  and  (b),  for  high  R/E  ratio  (R/E  = 
300).  (E/C  =  1.0,  N  =  25,  mesh  dimension  =  10^  x  10^  for  both  cases). 
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Figure  11 


Relative  concentrations  of  various  reaction  species:  (a),  for  high  evapora¬ 
tion/condensation  ratios  (E/C  =  10.0);  and  (b),  for  low  E/C  ratio  (E/C  = 
0.1).  (R/E  =  3,  N  =  10  for  both  cases). 


Figure  12  Simulated  microstructure  for  low  evaporation  to  condensation  ratio  showing 
unreacted  silicon  and  pore  distribution  for  case  given  in  Fig.  10  (a). 

Figure  13  Simulated  microstructures,  (a),  grain  size  and  pore  distributions  and  kinetic 
curves  of  relative  concentrations  of  reaction  species  for  high  number  nuclei  per 
particle  (N  =  10);  and  (b)  for  low  numbers  of  nuclei  per  particles  (N  =  1) 
(R/E  =  3;  E/C  =  1  for  both  cases). 

Figure  14  Simulated  microstructure  with  density  of  nuclei  N  =  7.3  (R/E  =  300,  E/C  = 
!)• 

Figure  15  Simulated  microstructure  for  high  reaction  evaporation  ratio  and  low  density 
of  nuclei:  (a),  final  microstructure,  and  (b),  same  as  (a)  with  overlay  of  initial 
particle  distribution.  Notice  the  large  pores  in  region  A  of  high  initial  packing 
density,  and  in  region  B  denuded  of  initial  particles. 
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Figure  1 


Representative  microstructure  of  reaction  bonded  silicon  nitride. 


Radial  Distribution  Function 


Figure  2 


Radial  distribution  function  of  coordinates  for  centers  of  spheres  obtained  from 
a  molecular  dynamics  simulation  of  the  structure  of  dense  random  packing  of 
spheres,  using  the  approach  of  Finney  [19]. 
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Figure  3 


Distribution  of  dense  randomly  packed  spheres  in  3-D  obtained  by  a  molecular 
dynamics  program,  subsequently  sampled  by  a  plane,  giving  a  2-D  starting 
structure  for  the  RBSN  reaction  simulation. 
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Figure  4  The  overlapping  of  the  2-D  arrangement  of  random  circles  by  a  hexagonal  mesh 

defining  the  cells  for  the  cellular  automata  reaction  simulation.  The  material 
inside  the  circles  represents  particles,  black  cells  are  Si^N^  nuclei. 


Real  Microstructure  Simulated  Microstmcture 


Figure  5  Pore  size  distributions  for  real  and  simulated  microstructures.  The  simulated 
structure  had  reaction  parameters  as  follows:  R/E  =  0.3,  E/C  =  1,  N  =  25 
(initial  density  of  nuclei). 
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Figure  6  Distributions  of  pore  shape  parameters  for  real,  (a),  and  simulated,  (b),  mi¬ 
crostructures  of  RBSN.  The  simulated  structure  had  reaction  parameters  as 
follows:  R/E  =  0.3,  E/C  =  1,  N  =  7.6. 


(a)  Simulated  history  of  completion  of  the  RBSN  reaction  as  a  function  of  time 
(x  is  the  total  transformed  fraction)  showing  Johnson-Mehl-Avrami  kinetics; 

(b)  experimental  thermo-gravimetric  results  of  Sheldon  and  Haggerty  [1]  for 
the  RBSN  reaction.  The  four  curves  are  for  four  separate  observations. 
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Figure  8  Simulated  microstructures  for  low  (R/E  =  0.3):  (a),  high  (R/E  =  300);  (b),  ra¬ 
tios  of  reaction/evaporation  rates,  and  (c),  showing  positions  of  initial  random 
packed  particles.  Packed  particles  were  equal  diameter  discs  packed  in  2-D  to 
illustrate  agglomeration  effects  and  the  tendency  of  the  high  R/E  simulation 
to  “heal”  packing  defects.  (E/C  =  1,  Nuclei  density,  N  =  10.) 
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Figure  8(b) 
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Figure  9  Simulated  RBSN  microstructure:  (a),  for  R/E  =  0.3,  E/C  =  1,  N  =  25 
(b),  for  R/E  =  300,  E/C  =  1,  N  =  25; 
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Figure  10  Relative  concentrations  of  various  reaction  species:  (a),  for  low  reac¬ 
tion/evaporation  ratio  (R/E  =  0.3);  and  (b),  for  high  R/E  ratio  (R/E  = 
300).  (E/C  =  1.0,  N  =  25,  mesh  dimension  =  10^  x  10^  for  both  cases). 
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Figure  11  Relative  concentrations  of  various  reaction  species;  (a),  for  high  evapora¬ 
tion/condensation  ratios  (E/C  =  10.0);  and  (b),  for  low  E/C  ratio  (E/C  = 
0.1).  (R/E  =  3,  N  =  10  for  both  cases). 


Figure  12  Simulated  microstructure  for  low  evaporation  to  condensation  ratio  showing 
unreacted  silicon  and  pore  distribution  for  case  given  in  Fig.  10  (a). 
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Figure  14  Simulated  microstructure  with  density  of  nuclei  N  =  7.3  (R/E  =  300,  E/C 
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Figure  15  Simulated  microstructure  for  high  reaction  evaporation  ratio  and  low  density 
of  nuclei:  (a),  .inal  microstructure,  and  (b),  same  as  (a)  with  overlay  of  initial 
particle  distribution.  Notice  the  large  pores  in  region  A  of  high  initial  packing 
density,  and  in  region  B  denuded  of  initial  particles. 
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SIMULATION  OF  CRACK  PROPAGATION  RESISTANCE 
IN  BRITTLE  MEDIA  OF  HIGH  POROSITY 

F.  G.  Haubensak,*  V.  V.  Bulatov,  and  A.S.  Argon** 
Massachusetts  Institute  of  Technology,  Cambridge,  MA  02139 

ABSTRACT 


The  crack  propagation  resistance  through  a  porous  or  microstructurally  heterogeneous 
brittle  solid  with  local  variability  in  strength  and  stiffness  has  been  simulated.  Specifically, 
the  simulation  probes  the  behavior  of  porous  brittle  materials  in  the  range  of  porosity  less 
than  cellular  materials  and  greater  than  microstructures  which  fall  under  the  category  of 
dilute  porosity.  The  simulation  plane  consists  of  a  triangular  network  of  points  interacting 
with  each  other  through  both  linear  central  force  springs  and  bond  angle  springs  incorpo¬ 
rating  an  appropriate  element  of  a  non-central  force  contribution.  Explicit  microstructural 
details  were  incorporated  into  the  model  and  the  simulation  was  first  carried  out  under  a 
condition  of  uniaxial  tensile  strain  in  order  to  investigate  the  mechanisms  of  sub-critical 
damage  evolution  leading  to  quasi-homogenous  fracture.  In  order  to  investigate  material 
strength  and  stiffness  variability  on  the  scale  of  a  representative  volume  element  for  coherent 
fracture  events  in  a  crack  tip  stress  gradient,  the  explicit  microstructural  results  were  incor¬ 
porated  into  a  simulation  with  boundary  conditions  characteristic  of  the  displacement  field 
of  an  infinite  Mode  I  crack.  To  impart  some  3-D  realism  to  the  primarily  2-D  simulations  a 
special  2-D  super  element  was  devised  which  incorporated  variability  information  as  might 
be  sampled  by  a  crack  front  in  3-D. 

*  present  address,  Department  of  Materials  Science,  Stanford  University,  Stanford  CA 

*•  author  for  correspondence. 
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For  a  given  porosity,  in  general,  only  small  differences  were  found  between  nominally 
diverse  microstructures  in  terms  of  their  tensile  toughness,  maximum  strength  and  elastic 
moduli.  The  strongest  dependence  of  the  overall  fracture  toughness  was  found  to  come 
from  the  average  porosity.  The  variability  in  local  element  strength  and  stiffness  on  the 
scale  of  the  porosity  produced  highly  tortuous  crack  paths  roughly  on  the  scale  of  the  chosen 
representative  volume  element.  The  tortuosity  of  the  crack  was  largest  where  local  variability 
of  strength  and  stiffness  was  uncorrelated.  Examples  of  microcrack  toughening  and  crack 
bridging  were  observed. 
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LIST  OF  SYMBOLS 


SYMBOL  DEFINITION* 

A  fracture  surface  area 

Ao  total  fracture  surface  area 

E  Young’s  modulus 

Eo  Young’s  modulus  for  fully  dense  material 

Gc  critical  energy  release  rate  of  porous  material 

Geo  critical  energy  release  rate  for  fully  dense  material 

K*  fracture  toughness  of  effective  2-D  super  element  for  3-D  application 

K{(T,k)  fracture  toughness  of  a  microstructure  having  element  strength  a  and 

stiffness  k 

N  number  of  nodes  in  simulation 

Af  population  size 

elastic  strain  energy  of  bond  extension 

elastic  strain  energy  associated  with  flexure  of  the  angle  between  bonds 
ij  and  ik 

Uo  equilibrium  bond  length 

b  empirical  exponential  pre-factor  relating  porosity  and  elastic  modulus 

bo  equilibrium  bond  angle  term 

Cijk  bond  angle  (watch  spring)  stiffness 
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coefRcient  of  variation  of  variable  x 


force  on  a  site  due  to  pairwise  and  3-body  interactions 

grain  boundary  bond  strength 

bulk  (matrix)  bond  strength 

energy  release  rate  of  RVE,  (=  wvu) 

population  average  energy  release  rates  of  RVE’s 

linear  spring  constant 

population  average  value  of  bond  spring  constant  k  in  large  assemblies 
peak  stiffness  of  a  microstructural  element 
average  porosity 

scalar  distance  between  sites  i  and  j 

size  of  the  RVE  chosen  to  derive  properties  of  effective  3-D  super  ele¬ 
ment 

standard  deviation  of  variable  x 
unit  vector  connecting  sites  i  and  j 

tensile  toughness,  fracture  work  per  unit  volume  of  a  microstuctural 

element  tested  in  tension 

area  per-site  of  the  triangular  spring  network 

area  per-bond  of  the  triangular  spring  network 

total  elastic  strain  energy 


e 

Cc 

V 


<Ti 

a 

a* 

o-c 

(^0 

X 


strain  in  bond  or  RVE 

rupture  strain  of  an  individual  matrix  bond 

terminal  fracture  strain  of  a  microstructure  element  or  RVE 

Poisson’s  ratio 

stress  tensor  at  a  microstructure  nodal  element  or  RVE 
stress  on  a  bond  connecting  nodes  i  and  j 
population  average  value  of  bond  strength 

peak  strength  of  an  effective  2-D  super  element  for  3-D  application 
peak  strength  of  individual  matrix  bond  in  the  porous  microstructures 
peak  strength  of  a  microstructure  element 
two  variable  correlation  parameter 


‘Properties  italicized  are  local  properties,  global  properties  are  not  italicized. 
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1  INTRODUCTION 


The  understanding  of  the  mechanical  behavior  of  heterogeneous,  and  more  specifically,  porous 
materials,  is  of  increasing  technological  interest.  Research  has  concentrated  mostly  on  materials 
with  dilute  porosity  [1],  and  on  highly  porous  cellular  structures  which  have  relative  densities  as 
low  as  0.3  [2].  For  each  of  these  extremes  of  porosity  there  are  mechanical  models  that  are  based 
on  assumptions  about  the  specific  forms  of  the  microstructure.  In  each  case,  extrapolations  can  be 
made  to  intermediate  porosities  with  some  limited  success.  However,  there  is  a  fundamental  lack  of 
understanding  of  behavior  of  materials  within  the  intermediate  range  of  porosities  (0.15  <  p  <  0.85), 
exhibited  by  reaction  bonded  ceramics. 

The  lack  of  understanding  of  the  fracture  behavior  of  materials  within  this  intermediate  range 
of  porosities  stems  from  the  fact  that  the  interaction  of  a  macro-scale  crack  with  heterogeneous 
microstructures  is  generally  beyond  the  scope  of  analytical  treatment.  For  a  planar  crack  in  the 
dilute  porosity  limit,  two  first  order  effects  reduce  the  fracture  toughness.  The  relationship,  Kc  = 
\/E(jc,  implies  that  the  toughness  is  related  directly  to  E,  the  elastic  modulus,  and  Gc,  the  critical 
energy  release  rate.  The  average  elastic  modulus  decreases  roughly  as  E  ~  Eoexp(— bp)  where 
Eo  is  the  Young’s  modulus  of  the  fully  dense  material,  p  is  the  average  porosity  and  b  is  found 
experimentally  to  be  typically  2-5  [3].  In  the  absence  of  other  effects  which  disperse  the  action 
away  from  the  parting  plane,  the  critical  energy  release  rate  Gc  for  a  planar  crack  is  proportional 
to  the  fracture  surface  area,  A  =  (1  —  p)Ao,  such  that  Kc  =  \/GcoEo(l  —  p)exp(— bp/2),  where  Geo 
is  the  critical  energy  release  rate  of  the  fuUy  dense  material. 

Any  crack  interaction  with  porosity  which  disperses  the  action  away  from  the  plane  of  fracture 
produces  non-planar  fracture  that  will  modify  this  relationship,  by  such  effects  as  microcracking, 
crack  bridging,  and  other  three  dimensional  shielding  effects.  As  porosity  increases  and  the  elastic 
interaction  based  on  pores  becomes  important,  the  approximation  of  isolated  pores  breaks  down. 
This  can  be  expected  to  occur  roughly  near  a  porosity  of  16%,  the  percolation  threshold  of  random 
porosity  in  three  dimensional  structures. 

Similarly,  cellular  structure  models  can  be  extended  with  some  success  into  the  intermediate 
porosity  regime  [2].  However,  the  fracture  of  cellular  materials,  especially  open  ceU  structures 
for  which  no  sharp  crack  tips  exist,  is  qualitatively  different  from  the  fracture  of  higher  density 
materials.  The  mechanism  of  crack  propagation  in  cellular  materials  in  tension  is  typically  a 
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succession  of  localized  breaks  of  individual  cell  ligaments  or  cell  walls.  Random  microstructures 
with  densities  higher  than  the  percolation  threshold  of  random  spheres  in  3-D  (16%),  no  longer 
require  structural  order  in  order  to  support  load.  Therefore,  this  limit  also  serves  as  an  approximate 
reference  point  for  the  limit  of  applicability  of  cellular  structure  fracture  models. 

The  material  of  particular  interest  in  the  present  study  is  reaction  bonded  silicon  nitride,  Si^N^ 
(RBSN),  with  a  porosity  roughly  in  the  range  of  0.25.  Figure  1  shows  a  typical  microstructure  of 
RBSN  interacting  with  a  macro-crack.  It  shows  relatively  equiaxed  pores  with  a  relatively  narrow 
pore  size  distribution  in  the  submicron  range.  In  addition  to  these  widespread  features,  RBSN  also 
contains  a  certain  small  concentration  of  unreacted  5i,  and  in  some  cases  larger  pores,  possibly 
resulting  from  initial  packing  defects  among  Si  particles.  Parallel  experiments  with  RBSN  revealed 
that  cracks  primarily  propagate  along  grain  boundaries,  bridging  the  pores,  as  is  shown  in  Fig.  1. 

In  a  related  companion  computer  simulation  [4],  we  have  presented  a  reaction  kinetics  model 
for  the  production  of  porous  SizN^  to  establish  the  important  reaction  parameters  that  give  rise 
to  microstructures  similar  to  the  real  microstructures  presented  in  Fig.  1,  and  to  identify  the  most 
important  factors  of  the  reaction  which  control  the  microstructures.  In  the  present  simulation 
the  fracture  processes  in  such  computer  generated  microstructures  were  explored  with  the  goal  of 
determining  the  cause-effect  and  mechanistic  relationship  between  features  of  the  microstructures 
and  their  toughness  which  can  be  translated  to  the  real  RBSN  system,  and  can  perhaps  guide  the 
further  process  developments  for  this  material. 

The  approach  taken  here  is  to  examine  the  fracture  response  of  realistic  porous  microstructures 
of  RBSN  in  as  much  detail  as  possible.  The  microstructural  details  are  incorporated  into  numerical 
simulations  by  mapping  them  onto  appropriate  degrees  of  freedom  of  the  spring  network  model,  as 
described  in  Section  2.  Also  discussed  in  Section  2  is  the  implementation  of  numerical  solutions 
and  other  pertinent  simulation  issues.  In  Section  3  we  begin  with  simulations  of  the  fracture 
behavior  of  brittle  porous  microstructures  of  RBSN  modeled  in  detail.  The  results  discussed  in 
Sections  3.1  and  3.2  show  that:  (1)  reliable  estimates  of  the  fracture  toughness  of  RBSN  from 
the  detailed  simulations  of  the  crack  propagation  within  the  spring  network  model  in  a  sample  of 
reasonable  size  are  hardly  feasible  because  of  the  high  computational  cost;  (2)  evolution  of  damage 
in  the  porous  microstructures  involves  sequences  of  fracture  cascades  leading  to  bridging  of  the 
inter-pore  ligaments.  While  the  first  observation  is  rather  discouraging,  the  second  implies  that 
substantial  computational  effort  can  be  saved  by  using  a  less  detailed  description,  where  the  porosity 
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is  incorporated  not  explicitly  but  in  the  form  of  variations  in  local  strength  and  stiffness  on  the 
scale  of,  pore  clusters  in  a  homogenized  model.  The  results  of  simulations  of  crack  propagation  in 
such  homogenized  materials  with  prescribed  variability  in  local  properties  are  presented  in  Section 
3.3.  The  first  two  sub-sections,  3.3.1  and  3.3.2,  describe  fracture  patterns  and  qucditative  trends 
in  fracture  toughness  observed  in  simulations  of  spring  network  models  with  various  idealized 
distributions  of  local  strength  and  stiffness  ranging  from  very  large  variations  to  no  variations 
at  all  in  one  or  both  local  properties,  and  from  no  correlation  to  perfect  correlation  of  them.  These 
results  demonstrate  the  range  of  changes  in  the  fracture  properties  as  functions  of  the  local  property 
variations.  To  incorporate  more  realistic  variations  and  co- variations  of  local  properties  of  porous 
RBSN  microstructures  in  the  homogenized  description,  stiffnesses  and  strengths  were  sampled  on 
large  statistical  ensembles  of  super  elements  of  the  size  10  x  10  corresponding  to  the  typical  effective 
pore  clusters  (Sub-section  3.3.3).  Then,  this  information  was  incorporated  in  crack  propagation 
simulations  of  the  homogenized  model,  now  representing  the  porous  RBSNs.  The  results  discussed 
in  Sub-section  3.3.4  show  that  local  variations  in  stiffness  and  strength  in  RBSNs  are  strongly 
correlated  in  all  cases.  The  origin  of  this  correlation  is  in  the  porosity,  which  largely  defines  the 
overall  fracture  properties  of  the  porous  microstructures.  Although  interesting  trends  are  observed 
and  discussed  in  Section  4,  the  nature  of  the  RBSN  formation  reaction,  limits  the  possibilities  of 
designing  tougher  materials  by  changing  the  reaction  conditions. 
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2  MODELING  FRACTURE  IN  A  BRITTLE  POROUS  MATERIAL 


2.1  Outline  of  the  Numerical  Model 

The  model  for  the  elastic  response  of  a  solid  in  two  dimensions  that  has  been  used  in  the 
simulations  described  here  is  based  on  the  properties  of  a  triangular  network  of  springs,  and  is 
similar  to  the  models  used  first  by  Kirkwood  [5]  and  by  Keating  [6]  and  more  recently  by  Srolovitz 
and  co-workers  [7].  The  network  is  composed  of  a  two  dimensional  triangular  lattice  of  points 
connected  by  bonds.  These  bonds  combine  two  different  forms  of  interactions:  first  a  two-point 
interaction  in  the  form  of  linear  harmonic  springs,  and  a  three-point  interaction  in  the  form  of 
angular  springs  referred  to  as  “watch”  springs  [7],  associated  with  the  angles  formed  by  the  adjacent 
linear  springs  originating  from  the  network  nodes,  as  shown  in  Fig.  2.  This  triangular  network  bond 
model  of  interconnected  extensional  springs  and  their  angular  flexing  interactions,  as  used  here,  is 
adequate  for  obtaining  numerical  solutions  for  various  problems  in  plane  isotropic  elasticity.  In  the 
form  the  procedure  is  used  by  us,  it  is  well  suited  to  provide  stress  and  associated  displacement 
field  solutions  for  systems  incorporating  dispersed  heterogeneities. 

In  brittle  materials  with  quasi-homogeneously  dispersed  porosity  the  list  of  possible  fracture 
toughening  mechanisms  that  can  be  present  include  microcracking  and  other  crack  tip  shielding 
mechanisms,  such  as:  crack  deflection,  crack  branching,  and  crack  bridging  analogous  to  fiber 
reinforcement  by  locally  tough  elements.  These  mechanisms  do  not  involve  plasticity.  Therefore, 
in  the  present  model  we  exclude  any  dissipative  and  non-linear  effects  such  as  stress  induced  shear 
transformations,  fiber  puU  out,  and  large  strain  deformation,  and  consider  only  local  variations  in 
strength  and  stiffness  of  the  representative  components,  as  is  the  case  with  RSBN  for  which  the 
model  is  principally  intended. 

The  energy  $  of  the  2-D  triangular  network  is  made  up  of  a  sum  of  bond  stretching  and  angle 
flexing  forms: 

e  (1) 

•  j  i  j.(fc=j+l) 

where  r,y  is  the  distance  between  node  (or  site)  “z”  and  any  one  of  its  nearest  neighbor  nodes 
“j”,  C/®  represents  the  energy  of  bond  extension,  and  ff"  represents  the  energy  associated  with  the 
flexure  of  the  angle  between  bonds  “zj”  and  “ik”.  The  factor  1/2  in  the  first  term  accounts  for  the 
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double  counting  over  all  network  nodal  points  N,  in  the  entire  field.  The  bond  extension  and  angle 
flexing  components  of  the  elastic  strain  energy  can  be  given  specifically  as: 


(2) 


^ijk  ~ 


(3) 


where  Ujj  and  Uik  are  unit  vectors  pointing  from  the  node  “i”  to  its  nearest  neighbors  “j”  and  “A:”, 
(•)  represents  a  scalar  product,  Oq  and  bo  are  the  equilibrium  bond  length  and  the  cosine  of  the 
bond  angle  in  the  unstressed  network,  and  k  and  c  are  the  extensional  and  watch  spring  stiffness 
constants  of  the  stretch  and  angle  flexing  interactions  of  the  network  nodes. 

To  relate  the  spring  network  response  to  the  actual  material  elastic  properties  we  note  that  the 
isotropic  elastic  properties  of  a  homogeneous  network  can  be  given  by  its  Young’s  modulus  E  and 
Poisson’s  ratio  u,  and  that  these  can  be  expressed  as  functions  of  the  two  stiffness  constants  k  and 
c  [7].  Thus,  the  effective  local  isotropic  elastic  constants  E  and  i/  for  a  uniform  spring  network 
model  are  given  by: 


E  =  Zka 


22kal  +  7c 
6kal  +  7c 


u  = 


2kal  —  7c 
Qkal  +  7c' 


(4) 


(5) 


Therefore,  by  a  proper  choice  of  k  and  c  the  network  represents  the  2-D  response  of  a  linear  elastic 
isotropic  solid  in  the  full  range  of  Poisson’s  ratios  from  -1  to  1/3.  The  hexagonal  symmetry  assures 
isotropy  [8]. 

In  order  to  reach  mechanical  equilibrium  of  the  spring  network  model  corresponding  to  a  min¬ 
imum  of  the  energy  function  $,  the  conjugate  gradient  method  was  used,  which  is  known  to  be 
numerically  efficient  and  robust  [9].  This  technique  allows  to  systematically  adjust  the  coordinates 
of  the  lattice  cites  by  iteration,  so  that  in  the  end  of  the  iterative  procedure  the  net  force  on  each 
site  is  zero  within  a  pre-set  tolerance  and  the  minimum  of  $  is  reached.  Such  states  of  mechanical 
equilibrium  are  conveniently  characterized  by  a  local  stress  tensor  definable  for  each  lattice  cite 
(or  bond)  of  the  spring  network.  For  operational  purposes  we  have  used  the  well-known  atomic 
stress  tensor  introduced  first  by  Born  and  Huang  [10],  properly  modified  to  account  for:  (1)  the 
non-central  forces  due  to  the  last  term  in  equation  1,  and  (2)  the  periodic  boundary  conditions 
employed  in  some  of  our  simulations,  as  explained  in  detail  in  Appendix  I.  In  that  appendix  we 
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also  demonstrate  that  the  conjugate  gradient  energy  minimization  technique  together  with  the 
generalized  definition  of  the  local  stress  constitute  a  fully  versatile  method  for  solving  boundary 
value  problems  for  strain/stress  distributions  in  isotropic  plane  elasticity.  It  is  shown  that  the 
stress  fields  obtained  for  a  simple  test  case  using  both  the  conjugate  gradients  method  with  the 
spring  network  and  the  finite  elements  method  for  the  corresponding  continuum  approach,  closely 
match  each  other.  This  correspondence  and  the  ease  with  which  the  variability  of  micro-mechanical 
properties  of  RBSN  is  mapped  onto  the  spring  network  make  it  a  natural  choice  for  the  specific 
purposes  of  the  present  simulation. 

2.2  Method  of  Implementation  of  Solutions 

In  order  to  probe  the  elastic  properties,  the  tensile  strength  and  the  crack  propagation  resistance 
of  a  microstructure,  appropriate  local  elastic  and  strength  properties  were  assigned  to  the  spring 
network  that  overlays  the  field  of  heterogeneities  such  as  pores,  grains,  inclusions,  and  interfaces. 
The  microstructure  volume  element  was  then  subjected  to  a  set  of  uniaxial  border  displacements, 
and  was  equilibrated  internally  using  the  conjugate  gradient  method  until  the  overall  elastic  strain 
energy  function  #  was  minimized  for  a  unit  level  of  border  displacement  which  served  to  scale  the 
internal  displacement  field  distribution.  Forces  in  all  extensional  bonds  were  then  computed  and 
the  bond  subjected  to  the  highest  ratio  of  tensile  force  to  bond  strength  was  identified.  This  ratio 
was  then  used  to  re-scale  the  border  displacements  associated  with  the  probing  “unit”  imposed 
strain.  The  stiffness  of  this  bond  was  set  to  zero  to  account  for  its  fracture.  The  mesh  was  then 
unloaded,  reloaded,  and  equilibrated  to  identify  the  next  weakest  bond  that  can  be  broken,  and 
so  on,  until  the  entire  cross  section  of  the  microstructure  volume  element  was  severed.  To  obtain 
stress-straining  curves,  the  overall  system  stress  on  the  microstructure  was  found  as  the  area  average 
of  the  node  level  stress  tensors  (see  Appendix). 

The  fracture  criterion  was  chosen  as  a  maximum  tensile  bond  force  as  the  best  approximation 
to  a  critical  local  tensile  stress  over  other  criteria,  such  as  critical  energy  density  or  bond  extension, 
since  brittle  solids  fracture  when  local  stresses  reach  the  intrinsic  material  strength. 
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2.3  Fracture  Behavior  of  Individual  Porous  Microstructures 


In  the  first  simulation  the  fracture  behavior  was  determined  of  individual  microstructures  of 
a  size  substantially  larger  than  the  representative  volume  elements  (RVE)  of  the  porous  material. 
Parenthetically,  we  note  here  that  the  size  of  an  appropriate  RVE  should  be  chosen  as  a  small 
multiple  of  interpore  spacings,  over  which  smoothed  behavior  is  obtained,  as  we  will  discuss  in 
more  detail  in  Section  2.4  below. 

The  microstructures  used  in  the  simulation  were  generated  primarily  by  a  Monte  Carlo  reac¬ 
tion  kinetics  simulation  referred  to  in  Section  1  above  [4].  In  addition  we  have,  performed  some 
simulations  on  microstructures  obtained  by  digitizing  actual  micrographs  of  typical  RBSN  samples. 

The  specific  larger  scale  microstructure  elements  used  in  the  simulation  were  mapped  onto 
the  two  dimensional  triangular  network  in  the  following  manner.  Solid  elastic  elements  of  the 
microstructure  were  differentiated  from  pores  such  that  solid  bonds  have  full  stiffness  and  bonds 
falling  into  pores  have  zero  stiffness.  The  fracture  properties  were  incorporated  by  assigning  strength 
to  bonds.  Here  intergranular  bonds  were  differentiated  from  the  matrix  bonds  by  assigning  them 
a  lower  strength,  e.  g.,  /j  =  0.75  /^,  where  the  latter  stands  for  the  strength  of  a  bulk  (matrix) 
bond.  This  value  of  boundary  bond  strength  was  chosen  as  an  estimate  of  the  strength  necessary  to 
develop  the  experimentally  observed  prevalence  of  intergranular  fracture  in  polycrystalline  RBSN,  in 
a  manner  similar  to  what  was  done  by  Wang  et  al.  [7].  Widespread  intergranular  fracture  has  been 
observed  in  high  purity  RBSN  manufactured  by  Haggerty  and  co-workers  [11]  whose  research  has 
been  closely  linked  with  the  present  simulations.  In  the  present  model  the  dimension  of  the  bond  was 
taken  to  be  a  fraction  of  the  grain  size  in  order  to  capture  some  of  the  internal  deformation  gradients 
due  to  the  pore  and  grain  size  distributions.  At  the  same  time,  the  specific  fracture  criterion 
of  a  critical  bond  strength  used  here  approximated  the  actual  intergranular  fracture  behavior  of 
RBSN  [11],  as  already  stated.  Once  the  spring  mesh  was  assembled,  and  the  strengths  and  stiffnesses 
of  all  bonds  were  assigned,  the  spring  network  of  the  sample  microstructure  was  subjected  to 
uniaxial  tension  by  prescribing  border  displacements.  Periodic  boundary  conditions  were  imposed 
on  the  test  microstructure  elements  to  approximate  essentially  infinite  bodies  with  no  free  surfaces. 
The  size  of  the  mesh  was  64  x  64  nodes,  approximating  to  a  nominal  dimension  of  2  microns. 
This  was  usually  enough  to  incorporate  10-20  pores  of  0.1  micron  dimension,  i.  e.  a  field  of  roughly 
25  RVEs  (see  Section  2.4). 
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2,4  Variability  of  Strength  and  Stiffness 


To  explore  more  general  behavior  with  greater  ease,  without  considering  specific  porous  mi¬ 
crostructures,  simulations  were  performed  on  continuous  material  in  which,  however,  the  char¬ 
acteristic  variabilities  of  the  porous  microstructure  were  locally  introduced.  This  approach  was 
convenient  to  probe  the  effect  of  heterogeneities  on  the  crack  propagation  resistance.  Thus,  for  this 
purpose  a  second  set  of  simulations  was  performed  in  the  gradient  field  of  a  Mode  I  crack.  In  this 
case,  each  node  represented  the  variable  properties  of  a  volume  element.  This  was  accomplished  by 
introducing  a  series  of  increasing  bond  strength  and  stiffness  variations  randomly  for  each  bond  to 
obtain  some  idealized  trends  resulting  from  local  variability.  In  a  third,  subsequent  simulation  the 
question  of  the  dimension  of  the  volume  element  was  addressed,  by  performing  a  series  of  fracture 
simulations  on  specific  microstructures  that  could  be  considered  as  representative  volume  elements 
(RVE)s. 

In  the  second  simulation  the  mesh  was  pre-cracked  and  was  placed  in  the  strain  gradient  field 
of  an  infinite  Mode  I  crack  which  approximates  a  realistic  crack  tip  environment  in  a  solid  body  of 
large  dimensions.  Beyond  the  border  of  a  core  region  surrounding  the  crack  tip  the  material  was 
considered  homogenized.  Inside  the  core  region  the  variations  due  to  individual  bond  properties 
were  considered  explicitly  as  they  interacted  with  the  crack.  Care  was  taken  to  minimize  the 
possible  adverse  effects  of  the  borders  on  the  crack  tip  by  keeping  the  size  of  the  core  region  larger 
than  any  crack  tip  process  zone.  The  local  fracture  criterion  for  crack  advance  was  taken  as  a 
maximum  principal  tensile  stress  approximated  by  a  critical  tensile  bond  force,  as  already  stated 
above. 

To  explore  the  effect  of  variability  of  bond  stiffness  and/or  strength,  the  bond  properties  were 
randomly  assigned  values  between  a  maximum  of  1.0,  and  a  minimum  value  determined  by  a 
variance  which  was  defined  by  a  coefficient  of  variation  (i.  e.  as  the  standard  deviation  divided  by 
the  mean).  The  density  of  the  distribution  of  the  bond  properties  was  assumed  to  be  constant 
over  the  range  of  values  defined  by  a  minimum  and  the  maximum  value  (1.0).  In  this  manner,  the 
coefficient  of  variation,  c{x),  in  the  variable  x  is  defined  as: 
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where  x{—  (xmax  +  2:min)/2)  is  the  mean  of  the  distribution  in  which  Xmax  =  1-  Thus,  a  variable 
uniformly  distributed  on  the  segment  (0,1)  has  a  maximum  coefficient  of  variation  of  0.577.  In 
the  simulation  each  parameter  of  strength  and  stiffness  was  varied  either  independently,  or  the  two 
parameters  were  varied  together  with  a  prescribed  correlation,  established  from  specific  simulations 
of  properties  of  actual  microstructure  volume  elements,  as  discussed  below  in  Section  2.5. 

2.5  Incorporating  Microstructural  Information  into  Stress  Gradient  Simulation 


To  make  the  choice  of  variations  in  strength  and  stiffness  of  the  bonds  and  their  co-variation 
less  arbitrary,  possible  combinations  of  microstructural  stiffness  and  strength  were  sampled.  This 
was  accomplished  by  testing  a  large  number  of  small  microstructure  volume  elements  as  potential 
representative  volume  elements  (RVEs),  extracted  from  a  very  large  field  of  a  porous  microstructure 
by  a  random  selection  process.  These  microstructure  VEs  were  tested  under  uniaxial  stress  to 
obtain  their  initial  elastic  stiffness  and  their  peak  stress  during  fracture.  This  is  described  in  detail 
in  Section  3.3.4  below. 

Moreover,  to  introduce  an  element  of  3-D  realism  where  porosity  in  a  real  structure  is  randomly 
distributed  into  the  third  dimension  rather  than  being  columnar,  as  implied  in  a  2-D  model,  a  special 
two  dimensional  “super  element”  was  constructed  by  averaging  the  properties  of  three  neighboring 
RVEs  in  a  row,  chosen  randomly  from  a  very  large  computer-generated  microstructure  as  depicted 
in  Fig.  3.  The  upper  bound  average  strength  and  stiffness  properties  of  such  rows  of  3  elements 
then  constituted  the  “super  element”  with  the  desired  3-D  properties,  to  represent  an  arbitrary 
strip  of  a  3-D  crack  front,  projected  into  the  2-D  plane.  The  effective  strength  a*  of  the  “super 
element”  was  obtained  based  on  the  argument  of  Rose  [12]  that  the  average  energy  release  rate  of 
a  strip  of  crack  front  with  heterogeneities  is  approximately  equal  to  the  line  average  of  the  local 
energy  release  rates,  g,  of  the  neighboring  crack  front  line  elements.  As,  i.  e.  : 


cr*  = 


(7) 


K*  ^  7E(ELig.As)/3As  ^  jEYA=i9i 

V  37r?-„  ’ 

where  Tu  is  the  size  of  the  microstructural  RVE  used  in  the  preliminary  simulations  to  gather 
statistical  data,  gi  is  the  energy  release  rate  (area  under  the  stress-strain  simulation  of  the  individual 
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RVEs  times  of  an  individual  RVE),  that  this  approximation  is  reasonable  is  demonstrated  in 
Fig.  4  where  a  more  exact  result  obtained  by  Bower  and  Ortiz  [13]  for  a  crack  pinned  at  two 
tough  heterogeneities  and  bowing  into  a  region  with  lower  toughness  for  five  different  ratios  of  the 
toughnesses  of  the  two  regions  through  which  this  crack  advances. 

The  simulation  which  incorporated  the  2-D  super  element  was  a  final  variant  of  the  third 
simulation  in  which  correlated  strength  and  stiffness  information  of  actual  RVEs  was  considered. 
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3  RESULTS 


3.1  The  Evolution  of  Damage  Under  Uniaxial  Tension 

The  damage  evolution  in  homogeneously  stressed  microstructure  volume  elements  is  characterized 
by  cascades  of  individual  bond  fracture  events  with  the  fracture  process  almost  always  initiating  at  a 
pore  surface  and  developing  into  a  succession  of  breaks  which  comprise  the  evolution  of  a  microcrack. 
The  nature  of  fracture  was  predominantly  intergranular,  with  75-100%  of  the  breaks  occurring  along 
the  lower  strength  grain  boundary  bonds,  as  it  was  arranged  to  happen  by  the  deliberate  choice  of 
the  lower  strength  for  the  boundary  bond,  to  parallel  the  experimental  observations.  The  cascade 
of  breaks  commonly  represents  the  sequential  breaking  of  elements  in  alignment  between  two  pores, 
which  is  generally  the  smallest  dimension  between  two  adjacent  pores.  This  was  established  from 
the  observed  fact  that  the  number  of  breaks  in  the  cascade  correlates  directly  with  the  number  of 
mesh  elements  that  make  up  a  ligament  that  fractures.  These  fracture  cascades  were  separated 
often  by  single  breaks  at  pore  surfaces  which,  however,  did  not  grow. 

Figure  5  shows  in  considerable  detail  the  response  of  a  stressed  microstructure  volume  element 
obtained  from  the  digitized  image  of  a  real  RBSN  sample.  The  figure  identifies  the  connection 
between  the  features  on  the  stress  strain  curve  with  events  on  the  stressed  microstructure. 

During  the  simulation  when  a  local  break  is  obtained,  such  as  at  the  point  near  A,  between 
the  two  central  pores,  the  system  registers  a  slight  decrease  in  elcistic  stiffness.  At  this  stage,  for 
operational  purposes  the  system  is  unloaded  and  subjected  to  the  same  unit  border  displacement  to 
seek  a  new  level  of  internal  equilibrium,  and  is  then  reloaded  until  the  second  break  is  registered  (at 
point  A).  The  decreased  stiffness  (increased  compliance)  of  the  system  manifests  itself  by  an  elastic 
loading  line  with  reduced  slope,  connecting  the  origin  to  point  A  (not  shown  in  the  stress  strain 
curve).  We  note  that  the  second  break  at  point  A  in  the  newly  re-equilibrated  microstructure 
occurs  at  a  system  stress  less  than  the  first  break.  This  was  a  frequent  occurrence.  After  this 
break  the  continued  loading  of  the  system  (following  the  usual  re-equilibration)  produces  a  third 
break  which  now  exceeds  the  stress  of  the  first  break.  All  other  dotted  lines  falling  inside  the 
outer  envelope  of  the  stress-strain  curve  represent  similar  virtual  re-equilibration  and  reloading 
processes  probing  the  system  in  the  computer.  In  a  real  system  loaded  by  border  displacements 
or  by  imposed  external  strain  increments,  there  would  be  no  retraction  of  strain,  but  the  system 
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response  would  be  given  by  vertical  stress  drops  as  the  outer  envelope  stress-strain  curve  indicates. 
In  this  case,  aU  individual  fracture  events  represented  by  crosses  fading  inside  the  stress  strain 
curve  represent  unstable  fracture  cascades,  as  those  identified  by  A-I  in  the  stress  strain  curve  and 
located  correspondingly  on  the  microstructure. 

During  the  evolution  of  damage,  leading  progressively  to  increasing  concentration  of  fractured 
ligaments,  the  overall  system  modulus  systematically  decreases,  shown  by  the  ever  decreasing  slopes 
of  the  lines  connecting  the  origin  to  the  new  fracture  events.  This  progressive  decrease  in  stiffness 
of  the  structure,  accumulating  increasing  levels  of  breaks,  is  summarized  in  Fig.  6. 

3.2  Uniaxial  Fracture  Behavior  of  a  Series  of  Microstructure  Volume  Elements 

A  series  of  porous  microstructures  which  were  produced  under  systematically  varied  condi¬ 
tions,  with  the  associated  reaction  kinetics  simulation  [4],  were  probed  in  tension  for  their  quasi- 
homogeneous  fracture  properties.  The  different  computer  generated  microstructures  shown  in  the 
various  parts  of  Fig.  7  were  deemed  to  be  realistic  in  appearance  and  within  quantitative  measures 
resembling  real  microstructures  such  as  the  one  shown  in  Fig.  1.  A  typical  control  microstructure  of 
25%  porosity  and  its  stress-strain  behavior  is  shown  in  Fig.  7a.  The  remaining  six  microstructures 
in  Figs.  7b  -  7g  represent  in  order;  a  case  of  a  higher  silicon  nitride  nucleus  density(rapid  nitrida- 
tion  and  early  particle  impingement)  (7b);  one  with  a  higher  porosity  at  31%  (7c);  and  one  with 
a  lower  porosity  at  20.5%  (7d);  a  case  of  a  significantly  anisotropic  grain  growth  condition  (7e); 
a  microstructure  for  a  case  of  a  low  reaction  to  evaporation  ratio,  in  the  reaction  kinetics  simula¬ 
tion  [4]  producing  a  controlled  variation  in  porosity  distribution  (7f);  and  finally  a  case  involving 
the  addition  of  a  population  of  small  diameter  (one  tenth  of  the  normal  size)  silicon  particles  (7g). 
This  last  case  was  undertaken  as  a  special  simulation  to  compare  with  the  experimental  observation 
of  a  measured  fracture  toughness  reduction  arising  from  the  addition  of  smaller  Si  particles,  in  the 
initial  “green”  material  and  was  thought  to  arise  possibly  from  different  pre-nitridation  particle 
packing  effects.  Details  of  the  reaction  simulations  used  to  obtain  microstructures  shown  in  Fig.  7 
are  discussed  by  us  elsewhere  [4]. 

The  corresponding  stress-strain  curves  resulting  from  the  simulations  are  also  shown  in  Figs.  7a  - 
7g  above  their  respective  microstructures.  In  aU  cases  the  stress-strain  curves  presented  are  obtained 
by  averaging  fracture  responses  of  four  similar  microstructures. 
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In  the  final  stages  of  the  simulations,  when  the  connectivity  of  the  microstructurers  is  substan¬ 
tially  reduced  due  to  damage  accumulation,  the  number  of  conjugate  gradient  iterations  required  to 
re-equilibrate  the  spring  network  becomes  increasingly  large.  To  reduce  the  computational  effort, 
the  final  “tail”  part  of  the  simulation  was  obtained  by  extrapolation.  For  that  a  certain  functional 
form  was  used  to  describe  the  “tail”  response,  which  was  obtained  by  an  empirical  fit  to  a  few 
typical  fracture  responses.  This  same  termination  form  was  then  applied  to  all  other  simulations. 

Perusal  of  the  composite  stress-strain  curves  of  these  simulations  shows  many  similarities  in  their 
shapes.  Apart  from  the  case  of  the  anisotropic  microstructures  which  shows  markedly  increased 
tensile  toughness,  all  other  cases  shows  few  quantitative  differences  apart  from  the  primary  effect 
of  changes  in  porosity. 

Figures  8a  -  8c  summarize  the  principal  results  of  these  simulations  on  specific  microstructure 
volume  elements,  showing  the  dependence  on  porosity  of:  the  Young’s  modulus,  E,  (8a);  the  peak 
strength,  CTc,  (8b);  and  the  fracture  work  per  unit  volume  w  (the  so-called  tensile  toughness)  as,  (8c). 
In  the  latter  the  beneficial  effect  of  the  anisotropic  grain  growth  stands  out.  In  all  other  cases  the 
porosity  level  is  clearly  the  most  important  structure  characterizing  parameter,  and  other  effects 
are  relatively  minor. 

In  the  simulations  the  tensile  toughness,  w,  was  obtained  as  the  area  under  the  individual 
stress-strain  curves,  i.  e.  : 

r^m 

w  =  I  ade  (8) 

Jo 

where  w  is  the  energy  per  unit  volume  required  to  completely  fracture  the  microstructure  to  a 
state  of  zero  load  support  under  uniaxial  strain,  and  is  the  strain  for  complete  fracture  of  the 
microstructure.  Typical  energy  densities  normalized  by  the  energy  to  fracture  a  perfect  elastic  solid 
of  bulk  material,  ac€cl2,  are  in  the  range  of  0. 1-0.2;  as  Fig.  8c  shows.  The  terms  (Tc  and  Cc  are  the 
strength  and  strain  to  fracture  a  fully  dense  matrix  in  uniaxial  tension. 

The  restriction  of  the  uniaxial  strain  condition  of  loading  can  produce  additional  load  support, 
especially  at  the  end  of  the  simulation,  due  to  the  absence  of  stress  relaxation  perpendicular  to 
the  loading  direction.  However,  the  magnitude  of  the  energy  contribution  under  the  tail  for  strains 
between  0.5  and  1.0  Cc,  the  critical  strain  to  failure  under  uniaxial  strain,  is  not  a  major  portion  of 
the  work  of  fracture  and  scales  with  the  total  energy,  tutaU  =  (0.3  -  0.4)rutotal- 
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3.3  Fracture  of  a  Fully  Dense  Reference  Solid  with  Prescribed  Variability  Mimicking 
the  Porous  SisN^ 

3.3.1  Uncorrelated  strength  and  stiffness  variation 

The  properties  of  porous  microstructures  in  the  uniaxial  tension  simulations  discussed  above 
revealed  that  because  of  the  porosity  and  the  variability  introduced  through  lower  strength  bonds 
for  grain  boundaries  in  comparison  with  the  strength  of  the  bulk  matrix,  damage  accumulation  was 
strongly  clustered.  This  effect,  which  was  discussed  in  detail  in  connection  with  Fig.  5  manifests 
itself  in  the  form  of  cracking  cascades  and  implies  that  for  purposes  of  simulation  of  crack  growth 
the  material  can  be  divided  up  into  smaller  volume  elements  as  discussed  briefly  in  Section  2.4 
above.  We  have  followed  this  approach  in  the  second  and  third  simulations  involving  crack  growth, 
as  outlined  in  Section  2.4. 

Thus  in  the  second  simulation  which  we  will  discuss  in  this  section  we  have  considered  the 
material  to  be  free  of  pores  in  which,  however,  the  expected  clustered  growth  of  a  crack  was  arranged 
to  result  from  material  heterogeneities  in  either  strength  or  stiffness,  as  might  result  from  porosity 
in  a  real  material,  and  where  the  variabilities  in  either  property  were  taken  as  uncorrelated.  In  the 
third  simulation  the  strength  and  stiffness  of  the  heterogeneities  were  correlated  on  the  basis  of  a 
separate  uniaxial  tensile  simulation  of  the  fracture  behavior  of  RVEs  of  a  size  roughly  corresponding 
to  the  cracking  cascades  that  were  observed  earlier  in  the  first  set  of  simulations  in  connection  with 
Figs.  5  and  7a  -  7g.  The  results  of  this  third  simulation  are  discussed  in  Section  4.2  below.  In 
both  the  second  and  third  simulations  we  have  utilized  the  same  triangular  spring  network  model 
discussed  above,  where  however  the  individual  “bonds”,  in  a  generalized  sense,  were  now  interpreted 
to  represent  the  heterogeneities  in  the  clustered  cracking  behavior,  noted  in  the  fracture  response  of 
the  actual  larger  microstructures.  The  propagation  of  Mode  I  cracks  through  such  “homogenized” 
material  with  prescribed  variability  was  then  investigated  in  a  manner  described  in  Section  2.5 
above,  leading  to  the  findings  described  below. 

The  degree  of  variation  of  both  the  bond  strength  and  stiffness  was  found  to  affect  directly  the 
extent  of  microcracking,  crack  deflection,  branching  and  bridging  behavior.  For  small  variability, 
the  crack  path  is  planar.  The  crack  samples,  on  the  whole,  the  population  average  of  element 
strength,  a,  and  the  average  fracture  toughness  is  that  of  a  homogenized  structure  with  the  average 
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properties,  Kc{a,k).  Here  Kc  =  ao\/Trao/2,  the  peak  fracture  toughness  of  the  heterogeneities,  has 
been  taken  as  the  unit  of  toughness  measure,  where  cTo  is  the  peak  strength  of  the  distribution  of 
material  elements  of  size  Oq  (the  appropriate  RVEs),  now  treated  as  a  distinct  “bond”.  As  the 
variation  of  local  properties  increases,  the  crack  path  becomes  more  tortuous,  and  the  number 
of  redundant  breaks  increases.  This  trend  is  shown  in  Figs.  9a  -  9f  for  increasing  coefficients  of 
variation  in  “bond”  strength  dc  of  elements  and  in  Figs.  10a  -  10c  for  increasing  coefficients  of 
variation  in  “bond”  stiffness  ko  of  elements. 

The  results  of  these  two  simulations  with  separate  strength  and  stiffness  variability  are  summa¬ 
rized  in  Figs.  11a  and  11b  respectively.  In  the  case  of  local  strength  variability  in  a  material  with 
uniform  local  stiffness,  the  toughening  is  a  direct  result  of  the  observed  increase  in  crack  complexity 
associated  with  the  very  substantial  redundant  microcracking,  which  could  not  be  present  in  the 
homogeneous  reference  structure.  This  is  clear  from  Fig.  11a  and  shows  that  the  beneficial  effect 
of  relative  improvement  in  toughening  reaches  a  maximum  with  a  coefficient  of  variation  of  around 
0.3.  In  the  case  of  the  largest  strength  variation  with  c(ct)  =  0.52,  the  toughening  effect  normalized 
with  the  toughness  of  the  homogeneous  material  was  1.26,  and  the  advancing  crack  splinters  into 
a  cloud  of  microcracks.  Figure  12  shows  the  resulting  R-curve  behavior  for  this  level  of  strength 
variability  arising  solely  from  the  tortuosity  of  the  crack  profile,  the  extreme  preponderance  of 
redundant  microcracking  and  the  associated  crack  deflection.  The  ratio  of  the  maximum  steady 
state  toughness  to  initial  toughness  is  observed  to  be  approximately  3  for  the  chosen  random  dis¬ 
tribution  of  element  strength.^  By  comparing  this  behavior  with  the  results  in  Fig.  9  we  note  that 
this  apparent  beneficial  behavior  is  at  the  expense  of  a  lowered  ultimate  (plateau)  toughness. 

In  the  case  of  stiffness  variation  alone,  where  the  material  has  uniform  strength,  the  effective 
toughening,  summarized  in  Fig.  11b,  must  result  solely  from  crack  tip  shielding  attributable  to  the 
random  elastic  inhomogeneities  which  then  manifest  themselves  in  crack  deflection  with  very  little 
associated  redundant  microcracking.  Since  the  material  is  of  uniform  strength,  the  toughening  due 
to  stiffness  variability  increases  monotonically.  Here  the  average  stiffness  of  the  broken  bonds  is 
equal  to  the  average  of  the  entire  population,  suggesting  no  favoring  of  high  or  low  stiffness  bonds 
for  fracture.  Hence,  any  elastic  shielding  is  due  to  a  favorable  stress  redistribution  arising  from  a 

^In  this  simulation,  the  first  60  bonds  were  broken  one  at  a  time  with  intervening  steps  of  equilibration, 
and  the  following  400  were  broken  in  sets  of  20  highest  stressed  bonds  to  expedite  calculations,  thus  explaining 
the  sparcity  of  simulation  points  for  increasing  numbers  of  breaks. 


20 


broader,  and  in  this  case  only  downward,  distribution  of  bond  stiffness.  Thus,  in  comparison  with  a 
homogeneous  material  where  all  bonds  have  the  same  average  stiffness,  the  largest  toughening  effect 
is  1.34,  achieved  with  the  largest  coefficient  of  variation  of  0.52  in  bond  stiffness.  We  note  that  the 
homogeneous  reference  material  also  maintains  constant  toughness  by  the  argument  that  a  more 
compliant  bond  absorbs  proportionally  more  energy,  yet  the  surrounding  material  has  a  reduced 
elastic  modulus  which  is  related  to  the  bond  stiffness.  This  demonstrates  that  the  toughness  of  a 
homogeneous  brittle  material  whose  energy  release  rate  is  inversely  proportioned  to  its  stiffness  is 
independent  of  the  elastic  properties,  for  a  stress  based  fracture  criterion. 

The  increase  in  toughness  due  to  material  variability  manifests  itself  directly  in  the  form  of 
crack  complexity  or  degree  of  deviation  from  planarity  of  the  crack  as  is  already  clear  from  the 
insets  in  Figs.  9  and  10.  This  is  demonstrated  more  directly  by  the  calculated  root  mean  square  of 
the  vertical  positions  of  the  breaks  with  respect  to  the  median  plane  of  the  crack  which  make  up  the 
overall  crack  profile,  and  naturally  includes  effects  of  crack  tortuosity  and  redundant  microcracking. 
Figure  13  shows  the  positive  correlation  of  the  crack  path  deviation  from  planar  form,  with  the 
relative  toughness  AKfK.  The  parameter  AK  is  defined  as  the  difference  between  the  average 
K  sampled  locally,  minus  the  reference  toughness  K  =  ^Je  x  g,  of  a  homogeneous  material  with 
average  properties.  Here  the  term  g  stands  for  the  energy  release  rate  of  one  severed  element  (bond) 
{l/2)(T€ao,  where  a  is  the  local  element  (bond)  strength,  e  its  failure  strain,  and  Co  is  the  element 
dimension.  We  note  that  the  dependence  of  the  normalized  toughness  increment  AKfK  on  crack 
surface  roughness  is  steeper  for  the  case  of  stiffness  variability  where  the  material  is  of  uniform 
strength  and  the  toughness  increase  must  be  a  direct  consequence  of  crack  tip  shielding  due  to 
the  elastic  inhomogeneities.  Higher  complexity  in  crack  configurations  and  crack  tortuosity  due  to 
microcracking  is  a  direct  consequence  of  the  strength  variability  as  already  noted  above. 

3.3.2  Effect  of  correlated  strength  and  stiffness  variation  on  crack  propagation 

We  note  that  while  the  abstract  analysis  of  the  second  simulation  discussed  in  Section  3.3.1  of 
the  effects  of  uncorrelated  variability  of  the  strength  and  stiffness  was  quite  instructive,  in  most 
real  materials  the  strength  and  stiffness  of  material  elements  wiU  be  more  or  less  correlated.  This  is 
especially  the  case  with  porous  materials  where  pores  reduce  both  local  strength  and  stiffness.  Thus, 
in  preparation  of  a  more  complete  third  simulation  discussed  in  Section  3.3.4  below,  a  variant  of  the 
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second  abstract  simulation  was  performed  where  the  strength  and  stiffness  of  individual  elements 
were  varied  together  in  perfect  covariation  and  then  independently,  with  a  coefficient  of  variation 
of  0.346  for  strength.  The  resulting  effects  on  the  crack  profile  are  shown  in  Fig.  14  where  now 
the  average  fracture  toughness  in  the  material  with  correlated  strength  and  stiffness  variation  has 
decreased  to  K  =  0.65  from  the  reference  case  of  strength  variation  alone  with  K  =  0.87,  measured 
relative  to  the  fracture  toughness  of  the  perfect  material. 

The  expected  value  of  toughness  for  a  cooperative  effect  between  strength  and  stiffness  vari¬ 
ability  is  A''expected/.fi^'(CT,fc)  =  K{a  -f-  da, ko)l K (a , ko)  x  K{ao,k  +  dk)lK{ao,k)  =  1.25  x  1.09  or 
1.37  ±  0.04,  when  normalized  with  the  homogeneous  equivalent  material  toughness,  where  ko  and 
(To  are  the  peak  element  stiffness  and  strength,  correspondingly.  The  actual  result  for  perfect  co¬ 
variation  of  strength  and  stiffness  is  a  crack  configuration  of  reduced  toughness,  as  compared  to 
the  homogenized  equivalent  material,  K(a  +  da,k  +  dk)lK(d,k)  =  1.04  ±  0.01.  The  uncertainty  is 
estimated  as  the  standard  deviation  of  the  mean,  6  =  where  s  is  the  population  standard 

deviation  and  M  is  the  population  size. 

The  conclusion  is  that  the  perfect  correlation  of  strength  and  stiffness  effectively  transforms  the 
fracture  criterion  into  a  critical  strain  criterion,  with  some  effects  of  elastic  shielding.  These  effects 
produce  weaker  microcracking  and  crack  diffuseness  and  lower  toughness  than  either  toughening 
effect  of  strength  or  stiffness  variation  alone. 

The  opposite  case  of  independent  strength  and  stiffness  variation  produces  increased  tortuosity, 
and  higher  toughness  with  much  higher  levels  of  redundant  microcracking  of  the  strength  variability. 
The  relative  toughness  is  then  K{a  -\-  da,k  dk)l K{a,k)  =  1.34  ±0.01,  which  is  an  effect  very  close 
in  magnitude  to  the  expected  effect.  Qualitatively,  the  meanderings  of  the  crack  configurations 
appear  similar  between  the  two  cases  of  correlated  and  uncorrelated  variations,  however  with  the 
latter  showing  a  much  higher  level  of  redundant  surrounding  microcracking. ^  This  is  confirmed 
by  a  nearly  equal  crack  profile  height  (r.m.s.  height  of  3.0  compared  to  2.6  element  dimensions  for 
the  uncorrelated  and  correlated  cases  respectively).  The  increased  toughness  for  the  uncorrelated 
material  is  due  to  two  effects.  The  elements  with  high  stiffness  and  low  strength  serve  as  potent 
microcracking  sources,  and  those  with  low  stiffness  and  high  strength  have  a  significantly  higher 
energy  absorbing  capacity. 


^This  explains  the  shorter  crack  for  the  same  number  of  breaks 
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3.3.3  Realistic  distribution  of  correlated  strength  and  stiffness  at  the  representative 
volume  element  level 

In  the  light  of  the  abstract  exercise  of  the  effects  of  the  correlated  and  uncorrelated  strength 
and  stiffness  variations  presented  in  Sections  3.3.1  and  3.3.2  above,  a  separate  simulation  was 
performed  to  obtain  more  realistic  information  on  covariation  of  strength  and  stiffness.  This  was 
derived  from  explicit  microstructural  simulations  of  the  properties  of  small  volume  elements  as 
specific  microstructure  units.  For  this  purpose  a  large  (1000  x  1000  points)  computer  generated 
microstructure  was  sectioned  into  three  families  of  increasing  cell  sizes:  with  10  X  10,  20  X  20  and 
40  X  40  points.  A  very  large  number  of  such  cells  in  each  family  of  three  sizes  was  then  tested  in 
tension  as  potential  candidates  for  RVEs  to  determine  their  distributions  of  porosity,  strength  and 
stiffness.  These  distributions  are  shown  in  Fig.  15.  The  results  of  the  10  X  10  point  cells  were  taken 
as  the  preferred  size  of  the  RVEs  since  this  came  close  to  the  extent  over  which  unstable  cracking 
cascades  had  been  observed  in  connection  with  the  first  simulation  discussed  in  Sections  3.1  and 
3.2.  In  this  manner  the  criterion  of  local  fracture  could  be  arranged  to  match  the  consequential 
material  response. 

Since  the  crack  propagation  resistance  of  various  microstructures  was  governed  directly  by  the 
local  variability  of  strength  and  stiffness,  and  their  covariation,  the  dependence  of  such  variability  on 
porosity  was  studied  in  lO'*  RVEs  of  10  X  10  size.  This  dependence  is  shown  in  Fig.  16  together  with 
the  specific  direct  correlation  of  the  strengths  and  stiffness  of  these  elements.  The  same  information 
on  variability  of  the  strength  and  stiffness  of  these  10^  RVEs  in  terms  of  their  coefficients  of  variation 
is  summarized  in  Table  I.  The  correlation  between  the  strength  and  stiffness  values  was  measured 
here  by  the  two  variable  correlation  parameter  x  as  follows: 

^  ^  [(o-(f0  -  a)(fc(fO  -  fe)] 

SaSk 

where  a  and  k  represent  the  strength  and  stiffness  of  the  RVEs,  r  is  the  location  of  the  element  in  the 
2-D  field  from  which  the  information  was  extracted  and  and  Sk  are  the  standard  deviations  of  the 
respective  distributions  of  strength  and  stiffness.  The  computed  values  of  x  were  found  to  be  in  the 
range  of  0.84  —  0.93  for  these  simulated  microstructures.  A  value  of  1.0  implies  perfect  covariation 
and  0  independent  populations.  The  high  value  of  x  implies  a  tendency  towards  the  critical  strain 
criterion  response  with  the  resulting  reduced  toughness.  In  addition,  the  variability  of  strength 
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and  stiffness  axe  seen  to  be  relatively  insensitive  to  the  details  of  the  different  microstructures  that 
were  generated  in  this  study.  The  coefficients  of  variation  of  these  RVEs  are  relatively  large  in 
comparison  with  those  considered  previously  in  Section  3.3.1.  They  range  from  0.59  to  0.68  and 
from  0.56  to  0.77  for  strength  and  stiffness  respectively.  The  strongest  dependence  of  variability 
was  seen  to  be  on  the  average  porosity;  as  porosity  increases,  the  variability  increases  for  both 
strength  and  stiffness,  as  shown  in  Fig.  17.  One  reason  for  this  is  that  higher  porosity  elements 
have  a  larger  number  of  possible  ways  for  this  porosity  to  be  distributed  (high  entropy),  within  the 
restrictions  of  the  reaction  kinetic  simulation  model  which  produced  these  microstructures.  This 
result  suggests  that  microstructures  with  a  higher  frequency  of  high  porosity  elements,  e.g.  high 
average  porosity  microstructures,  have  the  potential  for  higher  local  variability  and  consequential 
increased  crack  tortuosity,  and  associated  higher  energy  release  rate. 

In  addition  to  the  close  correlation  between  local  strength  and  stiffness,  the  connection  between 
local  stiffness  and  porosity  on  the  scale  of  the  RVEs  indicates  a  small  and  relatively  constant  effect 
of  the  typical  factors  other  than  average  porosity  which  govern  elastic  properties  of  porous  solids 
such  as  pore  size,  shape  and  spatial  distribution.  This  same  result  was  obtained  in  Section  3.1  for 
larger  microstructure  elements,  implying  that  over  the  range  of  scale  from  the  interpore  ligament 
to  a  scale  5-10  times  larger,  the  simulated  microstructures  exhibit  limited  sensitivity  of  elastic 
modulus  to  pore  shape  and  spatial  distribution,  and  strong  sensitivity  to  average  porosity,  as  the 
single  most  important  parameter. 

The  strength  and  local  porosity  are  seen  to  be  less  correlated  than  stiffness  and  porosity  (com¬ 
pare  Fig.  16(a)  and  Fig.  16(b)).  The  previously  mentioned  factors  for  variability  of  stiffness  asso¬ 
ciated  with  pore  morphology,  and  their  spatial  distribution  together  with  the  increased  variability 
generated  by  reduced  grain  boundary  strength  are  the  major  contributors  to  overall  variability  of 
properties. 

The  effect  of  grain  boundary  strength  is  seen  to  be  an  important  factor  in  the  de-coupling 
of  local  strength  and  stiffness.  This  is  a  consequence  of  a  decrease  in  the  covariation  parameter, 
x{o-,k),  from  0.951  to  0.854  produced  by  a  reduction  of  grain  boundary  strength  to  75%  of  fuU 
strength  of  matrix.  The  reduction  of  all  grain  boundary  strengths  in  such  a  manner  was  noted  to 
produce  the  high  proportion  of  intergranular  fractures  in  the  tensile  fracture  simulation  of  large 
microstructures  discussed  in  Section  3.1.  This  effect  also  serves  to  de-couple  the  local  strength  from 
local  porosity  as  evidenced  by  the  decrease  of  another  covariation  parameter  from  0.800  to 
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0.680.  The  fact  that  the  average  local  strength  is  reduced  from  0.287  to  0.267,  but  the  resulting 
toughness  measured  explicitly  by  propagating  a  crack  through  a  microstructure  remains  relatively 
unaffected,  demonstrates  the  effectiveness  of  the  variability  of  the  local  strength. 

3.3.4  Crack  propagation  with  correlated  specific  local  strength  and  stiffness  using  the  2-D 
super  elements 

To  evaluate  the  eflFects  of  specific  correlations  of  strength  and  stiffness  at  the  local  RVE  level  on 
crack  propagation  resistance,  and  to  incorporate  a  measure  of  three  dimensionality  to  the  microscale 
crack  propagation,  a  final  simulation  was  performed  utilizing  the  2-D  super  element  with  effective  3- 
D  properties  as  described  in  Section  2.5  above.  The  specific  correlations  of  strength  and  stiffness  of 
the  10  X  10  point  ceU  RVEs  discussed  in  Section  3.3.3  above  were  incorporated  into  the  properties  of 
the  “bonds”  in  a  random  manner.  Figure  18a  shows  the  effective  strength  and  stiffness  distributions 
of  the  2-D  super  cells  based  on  the  10  X  10  point  cell  RVEs,  averaged  3  at  a  time,  in  the  manner 
described  in  Section  2.5.  Clearly,  this  type  of  third  dimension  effect  has  removed  elements  of  zero 
strength  and  zero  stiffness,  compared  to  the  distributions  of  the  constituent  2-D  elements  shown  in 
the  upper  row  of  Fig.  15.  This  has  resulted  in  an  elevation  of  the  average  normalized  strength  to 
0.488  and  the  average  normalized  stiffness  to  0.365,  with  reductions  in  corresponding  coefficients  of 
variation  of  strength  and  stiffness  to  0.42  and  0.385  respectively,  from  the  values  given  in  Table  I. 
Figures  18b  and  18c  show  the  resulting  fracture  toughness  sampled  by  the  growing  crack  and  the 
tortuosity  of  the  crack  plane.  While  the  crack  plane  exhibits  significant  roughness  and  some  bridging 
behavior  there  is  minimum  redundant  microcracking  in  the  vicinity  of  the  crack  plane.  The  fracture 
toughness  normalized  with  that  of  the  homogeneous  equivalent  control  microstructure  (described  in 
Section  3.2)  is  now  1.15  ±0.02,  indicating  a  modest  effect  on  toughness,  compared  to  the  maximum 
possible  cooperative  toughening  of  ±  da^k  ±  dk)/ K{a\k)  =  (1.34  x  1.26)  =  1.69. 

In  summary,  the  major  result  of  the  simulation  of  crack  propagation  through  microstructures 
with  element  properties  obtained  by  sampling  RVE  scale  sections  of  simulated  microstructures  is 
that  the  measured  toughness  is  relatively  unaffected  by  variables  associated  with  pore  size  and 
shape  distributions,  but  most  strongly  determined  by  the  average  porosity.  The  simulated  fracture 
paths  were  of  a  similar  planar  character  with  minimal  microcracking  and  some  limited  bridging 
ligaments.  The  limits  of  pore  size  and  spatial  distribution  imposed  by  the  reaction  kinetics  model 


25 


used  to  generate  these  microstructures  serve  to  place  bounds  on  the  strength  and  stiffness  variations 
available  for  optimal  microstructure  design.  One  case  for  which  some  insight  can  be  derived  is  the 
high  porosity  microstructure  with  associated  relatively  higher  strength  and  stiffness  variabihty. 
This  case  demonstrates  the  competition  between  the  increase  in  toughness  due  to  higher  variability 
in  strength  and  stiffness  against  the  reduction  in  toughness  due  to  the  relatively  strong  dependence 
of  the  elastic  properties  on  average  porosity. 
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4  DISCUSSION  AND  CONCLUSIONS 


We  have  demonstrated  a  connection  between  pertinent  microstructural  features  and  crack  prop¬ 
agation  resistance  of  porous  heterogeneous  materials  by  simulating  the  fracture  process  in  computer 
generated  microstructures.  Such  models  demonstrate  that  by  associating  the  development  of  mi¬ 
crostructural  features  and  the  heterogeneities  in  local  stiffness  and  strength  with  processing  param¬ 
eters,  followed  by  probing  the  fracture  response  of  such  microstructures,  the  gap  can  be  bridged 
between  processing  and  mechanical  properties.  This  is  especially  relevant  for  fracture  toughness 
which  is  particularly  sensitive  to  microstructural  heterogeneities.  For  example;  the  variability  in 
local  strength  and  stiffness  affects  the  fracture  toughness  more  strongly  than  average  physical  prop¬ 
erties  such  as  the  elastic  modulus. 

The  two  dimensional  fracture  simulations  for  representative  microstructures  with  a  range  of 
20-30%  porosity  exhibit  a  stable  ligament  microcracking  behavior  which  results  in  the  evolution 
of  sub-critical  damage  leading  to  an  eventual  fracture  instability.  In  the  tensile  simulation  under 
displacement  control  no  instability  develops,  but  rather  characteristic  stress-strain  curves  with  long 
tails  of  decreasing  stress  are  obtained  for  material  degeneration  on  continued  extension.  The  tensile 
toughness  values  can  be  significantly  altered  in  comparison  with  those  of  a  perfect  homogeneous 
material.  This  tensile  toughness  ru  =  /  ade  ~  0.15(o’cec/2),  where  the  subscript  c  indicates  critical 
stress  and  strain  values  of  a  reference  perfect  material.  In  fact  an  otherwise  perfect  solid  with  a 
stress  amplifying  defect  will  have  less  energy  absorbing  capacity  than  the  porous  microstructures 
analyzed  here  if  the  stress  concentration  due  to  the  defect  is  more  than  approximately  3. 

The  frequency  and  effectiveness  of  beneficial  microcracking  are  strongly  determined  by  the 
blunting  ability  of  the  surrounding  pores  in  the  porous  microstructure.  For  three-dimensional 
microstructure,  the  possibility  exists  for  crack  growth  in  the  third  dimension,  and  the  blunting 
effect  may  be  diminished.  A  fuUy  three  dimensional  treatment  of  crack  interaction  with  porosity 
is  needed  to  fuUy  evaluate  this  effect.  An  attempt  in  the  introduction  of  a  2-D  super  element  with 
some  aspects  of  the  variability  of  the  material  in  the  third  dimension  has  provided  some  interesting 
trends  and  important  insight  indicating  the  mutually  beneficial  interaction  of  strong  and  weak 
segments  along  a  crack  front. 

Variability  of  elastic  stiffness  and  fracture  properties  on  the  scale  of  pore  clusters  give  rise 
to  increasing  crack  tortuosity,  microcracking,  bridging  and  associated  R-curve  behavior.  When 
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varied  independently,  the  local  stiffness  and  strength  each  produce  maximum  relative  toughness 
increases  of  34  and  31%  respectively.  Independent  assignment  of  strength  and  stiffness  (correla¬ 
tion  of  0)  to  local  material  elements  produces  a  tougher  structure  in  comparison  -with  proportional 
variation  of  strength  and  stiffness  (correlation  of  1).  This  latter  distribution  of  element  properties 
effectively  corresponds  to  a  critical  strain  fracture  criterion  and  a  toughening  effect  significantly 
less  than  either  of  the  two  toughening  effects,  determined  independently.  Effects  such  as  reduced 
grain  boundary  strength  (at  constant  bond  stiffness)  tend  to  decouple  the  correlation  of  strength 
and  stiffness  of  representative  volume  elements  and  results  in  increased  toughening  in  reference  to 
a  homogeneous  microstructure.  Samphng  the  properties  of  large  size  simulated  microstructures 
revealed  a  correlation  between  the  microstructural  element  strength  and  stiffness  of  0.84-0.93,  be¬ 
tween  the  two  extremes  of  0  and  1,  but  closer  to  the  less  potent  toughening  of  the  critical  strain 
criterion  behavior.  This  strong  correlation  is  a  direct  consequence  of  the  porous  microstructure 
where  pores  reduce  both  strength  and  stiffness.  In  other  situations  where  no  such  geometrical 
connection  exists,  less  correlation  between  local  strength  and  stiffness  should  be  attainable  which, 
with  increasing  variability,  has  significant  toughening  effects,  as  is  well  known  in  composites  [15j. 

The  question  of  microstructural  design  to  coax  the  crack  into  bridging,  crack  trapping  and 
blunting  configurations  is  key  in  order  to  improve  the  toughness.  We  have  shown  the  potency  of 
the  fracture  of  isolated  ligaments  to  increase  energy  absorption  in  a  representative  microstructure  of 
a  relatively  high  porosity  brittle  medium  and  quantified  the  role  of  local  material  variability  in  terms 
of  crack  plane  roughness  derived  by  effects  such  as  microcracking,  crack  deflection  and  bridging  for 
an  advancing  crack.  It  is  clear  that  crack  bridging  is  one  of  the  more  effective  mechanisms,  which 
is  predominant  in  highly  heterogeneous  structures  in  which  the  crack  tip  stress  singularity  is  stiU 
high.  This  is  the  mechanism  actually  observed  in  associated  experimental  studies  of  the  fracture 
of  porous  RBSN  [16]. 
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Table  I  Strength  and  Stiffness  Variabilities  of  RVEs  of  10  x  10  Size 


Simulation 

Description 

strength 

cr/o-o 

stiffness 

k/ko 

porosity 

P 

c(a) 

c(k) 

1 

b 

orrelatic 

a-p 

)n 

k  -  p 

K/Ko 

Control 

0.27 

0.32 

0.25 

0.61 

0.68 

0.85 

-0.68 

-0.84 

0.65 

(1) 

0.23 

0.32 

0.25 

0.64 

0.66 

0.93 

-0.77 

-0.84 

0.50 

(2) 

0.21 

0.24 

0.31 

0.68 

0.75 

0.81 

-0.60 

-0.81 

0.67 

(3) 

0.32 

0.43 

0.19 

0.60 

0.56 

0.93 

-0.79 

-0.87 

0.82 

(4) 

0.22 

0.30 

0.25 

0.59 

0.63 

0.93 

-0.75 

-0.81 

0.57 

(5) 

0.25 

0.29 

0.28 

0.66 

0.77 

0.84 

-0.65 

-0.83 

0.65 

(6) 

0.26 

0.32 

0.24 

0.62 

0.68 

0.85 

-0.67 

-0.85 

0.65 

Control 
(no  grain  bo 

0.29 

undaries) 

0.32 

0.25 

0.64 

0.68 

0.95 

-0.80 

-0.84 

0.63 
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APPENDIX.  THE  SPRING  NETWORK  MODEL  COMPARED  WITH 
THE  FINITE  ELEMENT  METHOD  IN  THE  SOLUTION  OF 
LINEAR  PLANAR  ELASTICITY  PROBLEMS 

In  order  to  examine  the  capabilities  of  the  spring  network  model  to  adequately  account 
for  non-uniform  stress  distributions  we  present  here  a  comparison  between  the  stress  field 
solution  with  the  spring  network  model  against  a  corresponding  solution  obtained  with  the 
finite  elements  method.  Since  the  focus  in  the  present  study  is  on  porous  solids,  an  infinite 
two-dimensional  periodic  array  of  circular  pores  is  chosen  as  the  test  case. 

For  the  sake  of  comparison  and  in  order  to  calculate  the  stress  distributions  a  local  stress 
tensor  is  introduced  either  on  the  basis  of  a  network  node  or  a  network  bond  in  the  spring 
network  model.  For  that  the  usual  definition  [10]  is  used  for  the  stress  tensor  expressed  in 
terms  of  the  “on-site”  forces  f,  as  follows: 

i^(fi(8)ri-|-r,<8)fi),  (A.l) 

where  a  is  the  stress  tensor,  0  is  the  total  volume  of  the  body,  r,-  is  the  position  vector  of 
site  i,  and  the  summation  is  performed  over  all  sites  within  the  body,  where  the  factor  1/2 
compensates  for  the  pair-wise  double  counting. 

As  was  demonstrated  by  Lutsko  [17],  equation  (A.l)  should  be  modified  for  the  case  of 
periodic  boundary  conditions  and  general  non-central  forces,  as: 

1  ^  d 

Oo-  =  -  ®  r.y  -f-  r,,-  <8)  {A.2) 

where  $  is  the  potential  energy  expressed  as  a  function  of  the  site  coordinates  and  r,j  =  Ti—Tj. 
By  substituting  equation  (A.l)  for  $  into  equation  (A. 2),  one  obtains  for  the  spring  network 
model  the  following  expression: 

i  j{i)  \  k(i,j)  } 
where  tensors  t,j  and  iijk  are  defined  as  ^ 

^All  other  quantities  in  equations  (A. 4)  and  (A.5)  were  defined  in  Section  2.1. 
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iij  =  kij{rij  -  ao)uij  (g)  (A.4) 

tijk  =  Cijkiuij  •  Uifc  -  bo){[uik  -  {Uijk  ■  Uiit)uij]  (8)  u,_,-  +  Uij  (g)  [u,-*  -  {uij  ■  Uik)uij]}/rij.  (A.5) 

In  equation  (A. 3)  summation  is  performed  over  all  sites  i  of  the  triangular  spring  network, 
their  six  nearest  neighbor  sites  j{i),  and  two  common  nearest  neighbor  sites  k{i,j)  of  sites  i 
and  j,  respectively. 

In  order  to  examine  the  stress  distribution  inside  the  body  it  is  necessary  to  partition  the 
total  stress  (equation  A. 3)  among  the  sites  (or  the  bonds)  requiring  that  the  volume  average 
of  local  stresses  be  equal  to  the  total  stress: 


ntr  =  Yl  '  E  (^-6) 

sites  bonds 

The  per-site  volume  is  defined  as  the  area  of  the  hexagon  surrounding  each  site. 
Accordingly,  the  “per-bond”  volume  is  defined  as  the  bond  share  of  the  “per-site” 

volume,  i.e.,  =  u;f‘^®/3.  It  is  natural,  then,  to  define  the  “per-site”  stress  as 


—site 


and  the  per-bond  stress  as 


1 


E  ( E  )  > 

AO  \  / 


(-4.7) 


=  5353  ( ‘.V  +  5  E  (t«i.  +  {.,  J  )  .  (4.8) 

b  \  k{i,j)  ) 

so  that  by  performing  summations  in  (A.6)  the  average  stress  (equation  A. 3)  is  obtained. 

Local  stress  computed  according  to  equation  (A. 8)  was  compared  against  the  correspond¬ 
ing  finite  elements  solution  for  the  hexagonal  periodic  array  of  circular  pores  subjected  to 
uniaxial  load  (Fig.  19).  Small  differences  in  the  component  (perpendicular  to  the  loading 
axis)  seen  near  the  pore  surfaces  are  due  to  the  discreteness  of  the  spring  network.  Other 
than  that,  the  agreement  is  almost  complete. 
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Figure  Captions 


Fig.  1  A  typical  microstructure  of  reaction  bonded  silicon  nitride  (RBSN), 

showing  the  propagation  of  a  crack  through  the  microstructure. 

Fig.  2  A3  node  point  configuration  of  the  extension  springs  and  angle  flexing 

“watch  springs”  in  the  mesh. 

Fig.  3  Method  of  determination  of  the  strength  and  stiffness  properties  of 

2-D  super  elements  by  line  averaging  the  properties  of  3  adjacent  rep¬ 
resentative  volume  elements  from  a  typical  field  of  simulated  porous 
microstructure,  as  explained  in  the  text. 

Fig.  4  A  comparison  between  the  stress  intensity  based  upon  the  weighted  line 

average  of  the  energy  release  rate,  g,  and  the  more  accurate  solution 
from  Bower  and  Ortiz,  [13],  for  crack  trapping.  Top  frame  depicts 
five  critical  crack  front  shapes  penetrating  into  a  low  toughness  region 
between  two  high  toughness  regions  for  the  five  specific  ratios  of  critical 
fracture  toughnesses  of  the  low  and  high  toughness  regions.  The  lower 
frame  shows  the  important  effect  of  crack  front  support  by  the  tough 
regions  for  the  five  special  solutions  of  toughness  ratio  of  the  low  to 
the  high  toughness.  In  the  range  of  ratios  from  0.4  to  1.0  the  effective 
toughness  is  well  described  by  the  simple  mean  value. 

Fig.  5  Fracture  response  of  a  uniaxially  stressed  microstructure  volume  el¬ 

ement  of  RBSN.  The  unstable  fracture  cascade  events  shown  in  the 
stress-strain  curve  are  associated  with  the  sites  where  such  fracture 
cascades  occur  in  the  top  figure  of  the  computer  generated  microstruc¬ 
ture. 

Fig.  6  The  monotonic  decline  of  the  normalized  Young’s  modulus  of  the 

stressed  microstructure  of  Fig.  5  with  increasing  damage  accumulation. 
Note  the  drops  associated  with  pore  linkage  and  subsequent  nucleation 
of  a  microcrack  elsewhere. 
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Fig.  7 


Fig.  8a 


Fig.  8b 


Fig.  8c 


Fig.  9 


Fig.  10 


The  microstructures  and  uniaxial  stress-strain  responses,  to  fracture, 
of  a  series  of  7  specific  cases:  (a)  a  control  sample  of  25%  porosity; 
(b)  a  case  of  high  initial  5'z3A^4  nucleus  density;  (c)  a  high  porosity 
microstructure  with  31%  porosity;  (d)  a  low  porosity  microstructure 
with  19%  porosity;  (e)  a  microstructure  with  anisotropic  grain  growth; 
(f)  a  case  with  negligible  product/reactant  impingement;  (g)  a  case  with 
4%  by  volume  of  smaller  Si  particles  in  the  green  material  population. 

Dependence  of  the  normalized  Young’s  moduli  on  porosity  of  the  7 
simulations  presented  in  Fig.  7.  The  key  to  the  identification  of  the 
information  is  as  follows:  (c)  control  microstructure,  Fig.  7a;  (1)  is 
Fig.  7b;  (2)  is  Fig.  7c;  (3)  is  Fig.  7d;  (4)  is  Fig.  7e;  (5)  is  Fig.  7f;  (6)  is 
Fig.  7g. 

Dependence  of  the  normalized  maximum  tensile  strengths  on  porosity 
of  the  7  simulations  presented  in  Fig.  7.  The  key  to  the  identification 
numbers  is  given  in  the  caption  of  Fig.  8a. 

Dependence  of  the  normalized  tensile  toughness  on  porosity  of  the  7 
simulations  presented  in  Fig.  7.  The  key  to  the  identification  numbers 
is  given  in  the  caption  of  Fig.  8a. 

Crack  propagation  experiments  through  fully  dense  material  with  vari¬ 
ability  in  bond  strength  but  constant  maximum  stiffness.  The  coeffi¬ 
cients  of  variation  in  strength  c(a)  increase  from  a  to  f.  Initial  positions 
of  the  pre-cracks  are  marked  by  short  vertical  lines. 

Crack  propagation  experiments  through  fully  dense  material  with  vari¬ 
ability  in  bond  stiffness  but  constant  maximum  strength.  The  coeffi¬ 
cients  of  variation  in  stiffness  c{k)  increase  from  a  to  e. 
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Fig.  11a 

Fig.  11b 

Fig.  12 

Fig.  13 

Fig.  14 

Fig.  15 

Fig.  16 


Summary  of  results  of  simulations  presented  in  Fig.  9  giving  the  de¬ 
pendence  of  average  fracture  toughness  on  the  coefficients  of  variation 
of  bond  strength. 

Summary  of  results  of  simulations  presented  in  Fig.  10,  giving  the  de¬ 
pendence  of  average  fracture  toughness  on  the  coefficient  of  variation 
of  bond  stiffness. 

R-curve  behavior  of  a  structure  with  a  high  coefficient  of  variation  in 
strength  of  c((Tc)  =  0.52. 

The  correlation  between  rms  crack  profile  height  and  the  relative  tough¬ 
ness  indicates  crack  deviation  from  planar  morphology  relates  to  mea¬ 
sured  toughness  due  to  such  effects  as  crack  deflection  and  microcrack¬ 
ing. 

The  effect  of  perfect  correlation  between  bond  stiffness  and  strength 
on  fracture  toughness,  for  a  coefficient  of  variation  of  0.346,in  both 
strength  and  stiffness,  compared  with  a  case  in  which  variability  in 
strength  and  stiffness  is  uncorrelated. 

Strength,  stiffness  and  porosity  distributions  obtained  by  probing  po¬ 
tential  representative  volume  elements  of  a  large  simulated  microstruc¬ 
ture  which  spans  10,000  volume  elements.  Upper  row  for  elements  of 
10  mesh  points  on  the  side,  middle  and  lower  rows  for  elements  20  and 
40  mesh  points  on  the  side,  respectively. 

Spread  of  normalized  stiffness  (a)  and,  strength  (b)  with  porosity  in 
10^  VEs  of  10  mesh  points  on  the  side,  tested  to  sample  realistic  distri¬ 
butions.  Correlation  of  normalized  strength  with  normalized  stiffness 
(c). 
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Fig.  17  Dependences  of  the  coefficients  of  variation  in  strength  and  stiffness  on 

porosity  for  the  simulations  presented  in  Fig.  16. 

Fig.  18  Effective  three  dimensional  strength  and  stiffness  distributions  embod¬ 

ied  in  the  2-D  super  element:  (a)  frequency  distributions;  (b)  evolution 
of  crack  propagation  resistance;  (c)  crack  plane  tortuosity. 

Fig.  19  Stress  field  solution  comparison  between  a  finite  element  generated  so¬ 

lution  and  the  discrete  spring  network  model  for  an  infinite  hexagonal 
array  of  pores  with  spacing  to  diameter  ratio  of  4  under  uniaxial  strain. 
Loading  direction  (xi)  is  vertical.  Heavy  black  lines  represent  the  con¬ 
tinuum  FEM  solution,  shading  and  thin  lines  represent  discrete  spring 
stress  solution. 
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A  typical  microstructure  of  reaction  bonded  silicon  nitride  (RBSN), 
showing  the  propagation  of  a  crack  through  the  microstructure. 


Fig.  2 


%  i 


A  3  node  point  configuration  of  the  extension  springs  and  angle  flexing 
“watch  springs”  in  the  mesh. 


Method  of  determination  of  the  strength  and  stiffness  properties  of 
2-D  super  elements  by  line  averaging  the  properties  of  3  adjacent  rep¬ 
resentative  volume  elements  from  a  typical  field  of  simulated  porous 
microstructure,  as  explained  m  the  text. 


CRACK 

PROPAGATION 
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Fig.  4  A  comparison  between  the  stress  intensity  based  upon  the  weighted  line 

average  of  the  energy  release  rate,  g,  and  the  more  accurate  solution 
from  Bower  and  Ortiz,  [13],  for  crack  trapping.  Top  frame  depicts 
five  critical  crack  front  shapes  penetrating  into  a  low  toughness  region 
between  two  high  toughness  regions  for  the  five  specific  ratios  of  critical 
fracture  toughnesses  of  the  low  and  high  toughness  regions.  The  lower 
frame  shows  the  important  effect  of  crack  front  support  by  the  tough 
regions  for  the  five  special  solutions  of  toughness  ratio  of  the  low  to 
the  high  toughness.  In  the  range  of  ratios  from  0.4  to  1.0  the  effective 
toughness  is  well  described  by  the  simple  mean  value. 
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Fig.  5  Fracture  response  of  a  uniaxially  stressed  microstructure  volume  el¬ 

ement  of  RBSN.  The  unstable  fracture  cascade  events  shown  in  the 
stress-strain  curve  are  associated  with  the  sites  where  such  fracture 
cascades  occur  in  the  top  figure  of  the  computer  generated  microstruc¬ 
ture. 


Normalized  Elastic  Modulus 


Fig.  6 


Number  of  breaks 


The  monotonic  decline  of  the  normalized  Young’s  modulus  of  the 
stressed  microstructure  of  Fig.  5  with  increasing  damage  accumulation. 
Note  the  drops  associated  with  pore  linkage  and  subsequent  nucleation 
of  a  microcrack  elsewhere. 
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Fig.  7  The  microstructures  and  uniaxial  stress-strain  responses,  to  fracture, 

of  a  series  of  7  specific  cases:  (a)  a  control  sample  of  25%  porosity; 
(b)  a  case  of  high  initial  Si3N4  nucleus  density;  (c)  a  high  porosity 
microstructure  with  31%  porosity;  (d)  a  low  porosity  microstructure 
with  19%  porosity;  (e)  a  microstructure  with  anisotropic  grain  growth; 
(f)  a  case  with  negligible  product/reactant  impingement;  (g)  a  case  with 
4%  by  volume  of  smaller  Si  particles  in  the  green  material  population. 
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Fig.  8a  Dependence  of  the  normalized  Young’s  moduli  on  porosity  of  the  7 

simulations  presented  in  Fig.  7.  The  key  to  the  identification  of  the 
information  is  as  follows:  (c)  control  microstructure,  Fig.  7a;  (1)  is 
Fig.  7b;  (2)  is  Fig.  7c;  (3)  is  Fig.  7d;  (4)  is  Fig.  7e;  (5)  is  Fig.  7f;  (6)  is 
Fig.  7g. 
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Fig.  8c 


Dependence  of  the  normalized  tensile  toughness  on  porosity  of  the  7 
simulations  presented  in  Fig.  7.  The  key  to  the  identification  numbers 
is  given  in  the  caption  of  Fig.  8a. 
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^  Crack  propagation  experiments  through  fully  dense  material  with  vari¬ 

ability  in  bond  strength  but  constant  maximum  stiffness.  The  coeffi¬ 
cients  of  variation  i.,  strength  c(e)  increase  from  a  to  f.  Initial  positions 
of  the  pre- cracks  are  marked  by  short  vertical  lines. 


Fig.  10  Crack  propagation  experiments  through  fully  dense  material  with  vari¬ 

ability  in  bond  stiffness  but  constant  maximum  strength.  The  coeffi¬ 
cients  of  variation  in  stiffness  c{k)  increase  from  a  to  e. 


Fig.  11a 


Summary  of  results  of  simulations  presented  in  Fig.  9  giving  the  de¬ 
pendence  of  average  fracture  toughness  on  the  coefficients  of  variation 
of  bond  strength. 
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Fig.  lib  Summary  of  results  of  simulations  presented  in  Fig.  10,  giving  the  de-  ^ 
pendence  of  average  fracture  toughness  on  the  coefficient  of  variation 
of  bond  stiffness. 
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Fig.  14  The  effect  of  perfect  correlation  between  bond  stiffness  and  strength 

on  fracture  toughness,  for  a  coefficient  of  variation  of  0.346, in  both 
strength  and  stiffness,  compared  with  a  case  in  which  variability  in 
strength  and  stiffness  is  uncorrelated. 
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Fig.  15  Strength,  stiffness  and  porosity  distributions  obtained  by  probing  po¬ 

tential  representative  v. ’lume  elements  of  a  large  simulated  microstruc¬ 
ture  which  spans  10,000  volume  elements.  Upper  row  for  elements  of 
10  mesh  points  on  the  side,  middle  and  lower  rows  for  elements  20  and 
40  mesh  points  on  the  side,  respectively. 


Spread  of  normalized  st'ffnets  (a)  and,  strength  (b)  with  porosity  in 
10“*  VEs  of  10  mesh  points  on  the  side,  tested  to  sample  realistic  distri¬ 
butions.  Correlation  of  normalized  strength  with  normalized  stiffness 

(c). 
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)endences  of  the  coefficients  of  variatic 
osity  for  the  simulations  presented  in 


)n  in  strength  and  stiffness  on 
Fig.  16. 
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Fig.  18  Effective  three  dimensional  strength  and  stiffness  distributions  embod¬ 

ied  in  the  2-D  super  element:  (a)  frequency  distributions;  (b)  evolution 
of  crack  propagation  resistance;  (c)  crack  plane  tortuosity. 


Stress  field  solution  comparison  between  a  finite  element  generated  so¬ 
lution  and  the  discrete  spring  network  model  for  an  infinite  hexagonal 
array  of  pores  with  spacing  to  diameter  ratio  of  4  under  uniaxial  strain. 
Loading  direction  (xi)  is  vertical.  Heavy  black  lines  represent  the  con¬ 
tinuum  FEM  solution,  shading  and  thin  lines  represent  discrete  spring 
stress  solution. 


